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Thermal Diffusion in Multicomponent Mixtures including H2O-CO,-Salts from 
Molecular Simulations 


Abbas Firoozabadi, Rice University 


Thermal diffusion is known to affect species segregation and separation in a wide range of 
problems including in the subsurface. This presentation will cover the formulations that will 
facilitate proper measurements, molecular simulations, and large-scale calculations in relation to 
mixture including subsurface CO2-H20-Salt systems. 


In laboratory measurement of thermal diffusion factor, the effect of thermodynamic factor is 
neglected. In both binary mixtures as well as in ternary and higher mixtures, thermodynamic factor 
may have a large effect. 


In ternary and higher mixtures, the Fickian diffusion coefficients may be required in evaluation of 
thermal diffusion factors. Molecular simulations can be used in estimation of thermal diffusion 
factors. In molecular simulations, a temperature gradient is established through non-equilibrium 
which has a number of subtleties. The implementation should be made carefully. 


CO>2 can be found in brine in the subsurface. CO2 permanent storage in the subsurface is the most 
promising method to prevent emission to the atmosphere. The long-term distribution of COz in the 
subsurface brine is of high interest. CO2 and brines may form complicated molecular structures. 
Some aspects of the problem in non-isothermal conditions will conclude the presentation. 
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Thermodiffusion of polymer solutions in mixed solvents 


W. Kohler', D. Sommermann!, M. Schraml! 


‘Physikalisches Institut, Universitat Bayreuth, Germany, werner.koehler@uni-bayreuth.de 


Introduction 


Diffusion and Thermodiffusion in polymer solutions have 
revealed a number of unexpected characteristics, such as the 
molar mass independence of the thermodiffusion coefficient 
(Schimpf et al. 1989), the speeding-up of diffusion with 
increasing polymer concentration in the semidilute regime 
(Adam et al. 1985, Rauch et al. 2003), and the insensitivity of 
the Soret coefficient to the glass transition at high polymer 
concentrations (Rauch et al. 2002). Little is known about 
thermodiffusion of ternary mixtures with polymers, e.g., a 
polymer dissolved in a binary solvent as represented by one 
of the DCMIX4 samples (Mialdun et al. 2020) and studied in 
a nonequilibrium fluctuations experiment (Garcia-Fernandez 
et al. 2019). 


Thermodiffusion in ternary mixtures has extensively been 
studied during the DCMIX microgravity campaign both in 
space and in the laboratory using, besides the 
thermogravitational column method, different two-color 
optical techniques. It has become clear that a major obstacle 
in all these experiments is the difficult inversion of the 
frequently ill-conditioned optical contrast factor matrix. In 
our contribution we explain how the inversion problem can 
be avoided in the case of strongly asymmetric ternary 
mixtures, such as polymers or colloids in mixed solvents. 
These systems typically show bimodal dynamics with well- 
separated diffusion eigenvalues. A natural choice, which 
needs to be justified by the results, is the assignment of the 
fast mode to the solvent and the slow mode to the polymer 
dynamics. This automatically fixes the directions of the 
diffusion eigenvectors in the space of the two independent 
composition variables and, thus, the diagonalization 
transformation of the diffusion matrix. 


In our contribution we analyse experiments within the 
framework usually employed for ternary mixtures. We show, 
how the ternary problem can be reduced to an effective binary 
one and we determine the diffusion and Soret coefficients for 
all components both for the asymptotic steady state and for 
the individual modes. These results are compared to data for 
binary reference systems, where possible. Using the 
transformation recently proposed by Ortiz de Zarate (Ortiz de 
Zarate 2020), transformation-invariant Soret coefficients are 
obtained. 


Experiment 


Experiments were performed on solutions of polystyrene (PS) 
of M, = 4880 g/mol dissolved in the mixed solvent toluene 
(Tol) and cyclohexane (cHex). The composition in mass 
fractions was 0.040/0.480/0.480 (PS/Tol/cHex). 
Measurements were carried out by means of 2-color optical 


beam deflection (2-OBD). Further details are described in 
Ref. (Sommermann et al. 2022). 


Results 


Let us take PS (ci) and Tol (cz) as the independent 
concentrations and cHex (c3) as the dependent one. Fig. 1 
shows a sketch of the concentration changes of all three 
components during the fast and the slow mode. 


fj ty 


Figure 1: sketch of the concentration changes during the fast and the 
slow mode for the case of r = c3/c2 = 0.5. Figure from Ref. 
(Sommermann et al. 2022). 


The identification of the concentration shifts associated with 
the two modes leads to the diffusion eigenvectors 


—— 


aa)? @* eden? \ 1 


where r = c3/c2 is the concentration ratio of the two solvents. 
The diffusion eigenvectors yield the diagonalization 
transformation matrix and, with the diffusion eigenvalues 


D, , the diffusion matrix 


j D, 0 
D = 7 : Vv = bp=D, r 
0 Ds tir D, 


Since the two diffusion eigenvalues are rather different, the 
off-diagonal element D2; is of similar magnitude as the 
diagonal ones. 


The theory has been worked out in detail in Ref. 
(Sommermann et al. 2022). Starting with the four amplitudes 
M;, and the two diffusion eigenvalues extracted from a 2-OBD 
experiment, it allows to compute not only the diffusion matrix 
but also the total Soret coefficients of the three components 
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for the asymptotic steady state separation and also the partial 
Soret coefficients for the two modes individually. 


The biggest advantage of our pseudo-binary approach is the 
avoidance of the inversion of the two-color contrast factor 
matrix and the comparability of the errors to the ones usually 
observed for binary systems. This huge increase of accuracy 
is shown in Fig. 2 for the two primed Soret coefficients of the 
independent components PS and Tol. The scattered points 
show results from Monte Carlo simulations, where all 
experimental quantities have been subjected to realistic 
statistical errors. 


error ternary 

error pseudo-binary 
pseudo-binary 
ternary 


1 1.5 2 2.5 
s',,/ (10° Kk’) 
Figure 2: Monte Carlo simulations of the effect of experimental 
error on the primed Soret coefficients of PS (1) and Tol (2) for the 
cases of the conventional ternary evaluation and the new pseudo- 
binary method. Figure from Ref. (Sommermann et al. 2022) 


All results can be found in Ref. (Sommermann et al. 2022) 
and will be given in the presentation. As an example, Tab. 1 
shows the transformation-invariant Soret coefficients 
according to Ortiz de Zarate (Ortiz de Zarate et al. 2020) of 
the two independent components PS and Tol together with the 
partial Soret coefficients for the two modes. The agreement is 
perfect. 


20°C 25°C 30°C 35°C] o 
46.3 44.0 41.9 39.6] 2.5 


units 


Sr1 10-2 kK! 
Sr.2|107? K7!| -2.36 -2.24 -2.15 -2.09]0.15 


See |10- °K | «2.46 -2.33 2.24 -2.17| 0.4 
Sst |10-3K—!| 2.46 2.33 2.24 2.17) 0.4 
S30" }10-3K-!) 47.5 45.2 43.0 40.6] 2.5 


Table 1: Transformation-invariant Soret coefficients (top two rows) 
and partial Soret coefficients for the fast and the slow mode. 


Finally, the viscosity-scaled thermophoretic mobility of the 
polymer, which can be calculated from the partial Soret 
coefficient of the slow mode and the corresponding diffusion 
eigenvalue, agrees with the universal value for nDr observed 
for binary polymer solutions (Stadelmaier et al. 2009). 


Conclusions 


The inversion of optical 2-color experiments for the 
investigation of diffusion and thermodiffusion in ternary 
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systems has proven notoriously difficult without additional a 
priori assumptions. For certain systems with well-separated 
diffusion eigenvalues, reasonable assumptions about the 
directions of the diffusion eigenvectors essentially eliminate 
these problems and increase the accuracy for ternary systems 
to the one known from binaries. A natural attribution of the 
two modes seems feasible for well-separated dynamics as 
observed for polymers or colloids in binary solvents at 
moderate concentrations. With increasing polymer or colloid 
concentration, the interparticle interaction between entities of 
the large species becomes more dominant and the attribution 
of the diffusion eigenvectors more difficult. So far, only the 
dilute to moderately semidilute case has been worked out. 
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Introduction 


In recent years, there has been a growing interest in mass and 
thermodiffusion in multicomponent mixtures, with a 
particular focus on ternary mixtures. The application of 
optical techniques has significantly improved the quality of 
mass and diffusion experiments, starting with binary liquid 
mixtures and expanding to ternary mixtures. Additionally, the 
resurgence of interest in thermodiffusion is attributed to the 
availability of convection-free experiments. 


The first microgravity experiment related to 
thermodiffusion was conducted in 1992 during the EUREKA 
mission, where the Soret coefficient was measured in 20 
binary organic mixtures and aqueous solutions (Van 
Vaerenbergh et al., 1995). While the experiment was only 
partly successful, experiments related to measuring Soret 
coefficients have been regularly conducted in orbit since then. 


IVIDIL ’ 
{ 2009-2010 \ 


DCMIX1 


2012 
THN-IBB-C12 


DCMIX2 
2014 
Tol-MeOH-Ch 


DCMIX3 
2016 
TEG-Wat-EtOH 


SODI 


DCMIX4 -2019 
‘ Tol-MeOH-Ch, C60-THN-Tol, 


Poly-Tol-Ch 


Figure 1: (inside blue frame) The outline of the Soret 
experiments carried out on the ISS in the SODI (Selected Opticals 
Diagnostic Instrument. (outside blue frame) Future experiment. 


Figure 1 illustrates experiments conducted in the SODI 
(Selectable Optical Diagnostic Instrument) on __ the 
International Space Station (ISS), equipped with a two- 
wavelength interferometer to measure Soret and diffusion 


coefficients in ternary mixtures. 

The first experiment, ITVIDIL (Influence Vibration on 
Diffusion in Liquids), demonstrated the feasibility of 
conducting experiments on the ISS and paved the way for 
ternary experiments. 


The DCMIX project is perhaps the most important 
collaboration in this area in recent years, in which teams from 
various countries have worked with ESA and Roscosmos to 
measure ternary liquid mixtures aboard the ISS. The DCMIX 
project also provided a valuable basis for the current GIANT 
FLUCTUATIONS project, which aims to investigate 
fluctuations in ternary systems subjected to temperature and 
Soret-induced concentration gradients under microgravity. 


The outcome of the DCMIX project is not only the 
measured coefficients; the amount of fundamental 
knowledge gained is much more significant. Aside from the 
measured coefficients, the DCMIX project resulted in 
numerous fundamental discoveries and breakthroughs during 
the experiment preparation and result processing stages. 
These findings are separated into two categories: generic 
(global) and experiment-dependent. The ‘generic’ results 
present a kind of scientific breakthroughs related to the Soret 
experiment in ternary mixtures, but independent of the 
particular DCMIX experiment. The ‘experiment-dependent’ 
results are new and interesting findings specific to the 
DCMIX mixture under consideration.. Due to space 
limitation only ‘generic’ results are shortly presented here. 


‘Generic’ results 


1. The first ever benchmark between ground and 
microgravity results as well as the first benchmark 
for ternary mixtures (Bou-Ali et al., 2015). 


2. Understanding the importance of the condition 
number for ternary and higher mixtures (Shevtsova 
et al., 2011). 


Optical methods, such as phase-shifting interferometry 
utilized in SODI, detect concentration changes in a sample 
using measurements of the refractive index (n). The 
conversion requires contrast factors, which are derivatives of 
n with concentration and temperature. In ternary mixtures, the 
2x2 matrix composed of the contrast factors can be ill- 
conditioned, limiting the mixture compositions that can be 
measured. This also applies to the thermogravitational 
column, where the derivatives of the refractive index and 
density are considered. 
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3. Revelation and insight into the asymmetric nature 
of error propagation in the Soret experiment. 


a) The error bar forms a very elongated ellipsoid instead of an 
isotropic cloud around the solution (Mialdun et al., 2018). 


b) Not all three Soret coefficients in a ternary mixture are 
obtained with the same accuracy. The individual uncertain- 
ties depend on the orientation of the error ellipsoid within the 
Gibbs triangle and its projections onto the three axes (Triller 
et al., 2019). 


4. New methodologies for processing large amounts of 
raw data from the Soret experiment, even partially 
corrupted by laser instabilities. 


a) Suggestion of an optimal multi-step path to process results 
with a smaller uncertainty (Mialdun et al., 2018). 


b) A robust data evaluation method that maps the SODI 
experiment to an equivalent OBD experiment and does not 
rely on phase stepping (Sommermann et al., 2019). 


5. Introduction of the Soret vector for component 
separation and the new relationship between binary 
and ternary coefficients on the Gibbs triangle. 


a) The Soret vector on the Gibbs triangle elucidates the 
essential features of Soret-driven separation in ternary 
mixtures as a whole (Mialdun et al., 2021). 


b) Sign changes along the binary boundaries and the existence 
of a singular point in the ternary Gibbs diagram, where all 
Soret coefficients vanish simultaneously (Schraml et al., 
2021). 


The Soret vector, as defined in Figure 2, offers several 
advantages: (i) to predict the Soret sign of a ternary mixture 
from knowledge of the Soret coefficients in binary 
subsystems; (ii) to control the consistency of measured 
coefficients; (iii) to determine the regions and components 
causing the greatest separation; and (iv) to identify the regions 
where Soret separation is inaccessible for optical techniques 
or gravitationally unstable. 


The peculiarity of the singular point in the ternary Gibbs 
diagram is that a mixture of this composition will remain 
homogeneous in the presence of a temperature gradient. 


Conclusions 


The DCMIX project has significantly increased knowledge 
about ternary mixtures. Microgravity experiments have 
played a key role in that. 


Microgravity experiments define crucial experimental 
cornerstones that are complemented by a large number of 
ground experiments. 


Microgravity experiments permitted a consolidation of strong 
international science teams. The tests on-board the ISS were 
an essential motivation. 
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Figure 2: Concept of the Soret vector in a ternary mixture 
wich is defined by its component as St'=Sr)' e1+Sr2'e2 where 
ej is the unit vector along the axis w;. The coordinates of the 
blue (red) part of the vector corresponding to the 
concentration change due to temperature decrease (increase). 
The two half-vectors, symmetrically located around the mean 
concentration, do visually show how the concentration 
evolution with temperature difference. 
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Abstract 


Recently, several European scientists have contributed to a 
joint effort led by the European Space Agency aimed at the 
identification of the priorities for future research in space in 
the field of Soft Matter and Biophysics (Vailati et al. 2021). 


In this work we present an overview some of the aspects that 
make the investigation of diffusion in liquid mixtures in 
Space still a challenging and promising field of research, of 
strategic relevance for space exploration, due to its 
applicative importance for the processing of fuels, food, and 
other materials needed for the sustainability of long-term 
space travels. 


In the past and present years, the investigation of diffusion in 
liquid mixtures has largely benefitted from the opportunity of 
performing experiments under reduced gravity conditions 
using sounding rockets and Low Earth Orbit platforms, such 
as the International Space Station (Braibanti et al. 2019). 
Until now, diffusive mixing has been mainly investigated at 
the macroscopic scale. Its investigation at the mesoscopic 
scale is becoming increasingly important for the 
understanding of mass transfer in confined systems, such as 
complex fluids, soft matter, biological systems, porous media, 
and microfluidic devices. Microgravity conditions provide 
the opportunity of investigating the effect of external fields 
on diffusive mixing in the absence of the convective flows 
induced by buoyancy on Earth (Vailati et al. 2011), with 
applicative relevance to the processing of liquids in space 
(Vailati et al. 2020). 

Even though microgravity experiments are still of primary 
relevance for the fundamental understanding of diffusion in 
complex fluids, the significance of these studies is shifting in 
the direction of investigating diffusion in multi-component 
mixtures under hypergravity and reduced gravity conditions, 
a topic of great relevance for long-term space exploration 
missions, where the ability of processing materials in the 
liquid state will represent a strategic habilitating factor 
(Vailati et al. 2023). 
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Introduction 


In the present work we investigate the influence of variations 
of the aspect ratio of the cell, , on the patterns arising in a 
slightly inclined cylindrical container heated from below for 
binary mixtures. In the presence of the shear mechanism 
induced by inclination, superimposed to the thermodiffusive 
effects, oscillatory superhighway convection (SHC) patterns 
coexist with the large-scale shear flow (LSF) (Figure 1). SHC 
patterns show a number of parallel thermal lanes, each 
containing aligned coherent structures that counter-propagate 
in adjacent lanes. LSF patterns generate a linear concentration 
profile across the cell. Both families of solutions are stable. 


LSF SHC 
Figure 1: Temperature and concentration contour plots for 
coexisting stable LSF and SHC states. 


The oscillatory SHC patterns arise in the Soret regime, where 
concentration gradients contribute significantly to the 
dynamics, and were first observed experimentally (Croccolo 
et al. 2013) and later obtained numerically (Alonso et al. 
2018). Some of the complex patterns obtained in a T=5 cell 
were analysed in detail in Alonso et al. 2022. 


We will present and analyze the spatio-temporal features of 
the patterns appearing in the Soret regime for values of the 
aspect ratio of [=5.2, 5.3, 5.4. When the value of T is 
changed, SHC-like patterns continue to arise, but their spatio- 
temporal properties differ considerably: symmetric and non- 
symmetric periodic orbits, modulated solutions with two and 
three frequencies, and chaotic solutions are some of the 
patterns that we obtain in wide regions of the parameter space. 


Equations and numerical tools 


We consider the Boussinesq binary-fluid convection in a 
cylindrical cell of height H and radius R, inclined an angle a 
with respect to the horizontal. The cylinder is heated from 
below, with a temperature difference between the lids equal 
to AT. The whole boundary is impermeable and non-slip, with 
fixed temperature at the lids, while the lateral wall is assumed 
to be thermally insulated. The governing equations of the 
problem reflect the incompressibility condition, the mass and 
heat conservation laws, and the Navier-Stokes equations in 
the Boussinesq approximation. 


The resulting system of nondimensional equations depends 
on the inclination angle a, the aspect ratio [=R/H, and four 
dimensionless parameters, the Rayleigh number Ra, the 
separation ratio S, the Prandtl o and Lewis numbers t. 
Equations and boundary conditions are equivariant under the 
group of symmetries that contains the transformations {J, Ri, 
R2z, R3}. 1 stands for the identity, R; is a reflection with respect 
to the middle longitudinal vertical plane, R2 is a point 
symmetry with respect to the center of the cylinder, and R; is 
the composition of the previous transformations. The 
nondimensionalized system of equations and boundary 
conditions is solved numerically using the algorithm 
described in Mercader et al. 2010. More details about the 
derivation of the equations, symmetries, and the numerical 
method can be found in some previous works in the same or 
related configurations (Alonso et al. 2018, 2022). 


To gain insight of the spatio-temporal properties of the 
patterns obtained by DNS we use a post-processing tool 
known as higher order dynamic mode decomposition 
(HODMD) (LeClainche & Vega 2017). This method 
decomposes the thermal and concentration fields as Fourier- 
like expansions that can be obtained from a set of spatially 
discretized snapshots, which are portraits of the system for 
equispaced values of time. The use of HODMD allows to 
discard transient dynamics and facilitates identifying the 
spatio-temporal nature of the solutions 


Main results 


We take as reference values for the simulations those of the 
experiment of Croccolo et al. 2013, where a isobutylbenzene— 
n-dodecane at 50% weight fraction binary mixture is used. 
The nondimensional parameters for this mixture are 
S=0.13, t=0.011, and o=16. 


We focus here in the results corresponding to a ['~5.3 aspect 
ratio cell. Unlike in the ['~5 cell, no stable periodic solutions 
are obtained in this case. The SHC-like patterns that arise are 
modulated quasi-periodic or chaotic solutions composed of 9 
convective lanes, stable solutions with different number of 
lanes are not found in the explored region of the parameter 
space. In contrast, stable 8- and 9-lane periodic solutions are 
observed for the [=5.2 cell. Modulated 9- and 10-lane 
solutions coexist in the [=5.4 cell. 


Figures 2 and 3 show the time series in nondimensional time 
units (1 n.t.u. corresponds to 20s) for the kinetic energy of a 
few selected solutions obtained in the range Ra=(1596,1620). 
The time dependence of these solutions differs substantially. 
Figure 4 shows the averaged Nusselt number (top) and the 
frequencies associated to the pattern obtained by FFT 
(bottom) as a function of Ra. 
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Figure 2: Kinetic energy as a function of time for modulated SHC 
states obtained for different values of Ra in a [=5.3 cell. 
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Figure 3: Kinetic energy as a function of time: hysteretic behaviour, 
different types of modulated SHC can coexist (Ra=1619). 
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Figure 4: Averaged Nusselt number (top) and frequencies of the 
solutions (bottom) as a function of the Rayleigh number. Only the 
frequencies of the modulations are included, the large frequency 
(f=0.24) of the basic SHC pattern is not shown. 


Solutions for Ra=1596,1600,1614,1620 are exemples of 
different types of quasiperiodic orbits. For all of the them, the 
modulation spans for a long time compared with the period of 
basic SHC solution. The period of the fast frequency 
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associated to the basic SHC pattern is 7;=4.17 n.t.u.=83 s. For 
Ra=1596 the periods of the modulation is around 7;,,=100 
n.t.u.=0.5h, and it raises as Ra increases. For the largest 
modulations the periods are around 7,,;=666 n.tu.=3.7 h 
(Ra=1620). Figure 5 shows the phase portrait in two 
symmetric points (x,y) and (x,-y) for two types of quasi- 
periodic SHC solutions: a R)-symmetric QP orbit (Ra=1600) 
and a non-symmetric QP orbit (Ra=1614). 


Solutions for Ra=1598,1612,1618 and Ra=1602-1608 are 
chaotic orbits, and solutions such as the one obtained for 
Ra=1619 (third plot in Figure 3) add a still larger modulation 
to the one already present (T2=1333 n.t.u.=7.4 h). The 
system exhibits hysteretic behaviour for certain ranges of Ra 
values. As an example, Figure 3 includes two different stable 
solutions obtained for the same value of Ra, Ra=1619. 


Figure 5: Phase portrait in two symmetric points (x,y) and (x,-y) for 
two types of quasi-periodic SHC solutions: a symmetric QP orbit 
(Ra=1600, left) and a non-symmetric QP orbit (Ra=1614, right). 


Conclusions 


We present and analyse the spatio-temporal features of 
relevant patterns arising in an inclined cylinder in the Soret 
regime for aspect ratios around 7=5. The reconstructions of 
periodic SHC solutions provided by HODMD suggest that 
they originate from a Hopf instability of the basic steady 
pattern of longitudinal rolls. Subsequent bifurcations give rise 
to several types of modulated solutions of slowly varying 
amplitudes, which can be quasiperiodic and chaotic orbits. 
Lengthy computations are required to establish precisely the 
spatiotemporal nature of the SHC-like patterns. 
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Introduction 


Supercritical carbon dioxide (COz) mixtures are essential for 
numerous environmental, power and chemical industry 
applications. Although thermophysical properties of CO2 and 
some of its binary mixtures are well-known, the knowledge 
of their transport properties is insufficient. In fact, reliable 
data on diffusion, and particularly, on thermodiffusion of CO2 
mixtures are scarce, especially in the extended critical region, 
where anomalous behavior of these properties has been 
observed experimentally and with molecular simulation 
(Guevara-Carrion et al. 2019, Chatwell et al. 2021). 


Besides experimental techniques, molecular modeling and 
simulation is a powerful tool for the prediction of dynamic 
properties of fluids, but also for the understanding of 
microscopic interactions that lead to the observed diffusion 
and thermodiffusion processes. In this work, we present the 
results of a molecular simulation study on diffusion and 
thermodiffusion of seven binary supercritical carbon dioxide 
(scCOz2) mixtures with hydrogen, methane, ethane, isobutane, 
benzene, toluene and naphthalene at concentrations near the 
infinite dilution limit. These mixtures were studied along the 
9, 10 and 12 MPa isobars in the temperature range between 
290 and 340 K. Because of pronounced density fluctuations 
in the extended critical region, molecular dynamics 
simulations are particularly challenging and require extensive 
sampling and large system sizes. Equilibrium molecular 
dynamics simulations and the Green-Kubo formalism were 
applied, employing simple, rigid, non-polarizable molecular 
models. 


Fick and Maxwell-Stefan diffusion coefficients, as well as the 
related thermodynamic factor were investigated along the 9 
MPa isobar for three solute mole fractions (0.5, 1.0 and 1.5 
mol%). The Fick diffusion coefficient at infinite dilution was 
obtained by extrapolation of the intra-diffusion coefficient of 


xjol 


the solute i to the infinite dilution limit, Di; =D, 


Diffusion 


It was found that slight variations of temperature and pressure 
in the extended critical region may cause the density of these 
mixtures to halve, the shear viscosity to double and diffusion 
coefficients to rise strongly. This anomalous behavior of both 
thermodynamic and transport properties is particularly 
pronounced near the so-called Widom line. In contrast to the 
strong composition dependence of the Fick diffusion 
coefficient in the extended critical region, the Maxwell-Stefan 
diffusion coefficient shows a rather weak dependence on the 
solute mole fraction, cf. Fig. 1. This can be explained by the 
strong changes of thermodynamic factor, which account for 
the thermodynamic non-ideality of a mixture and is, per 


definition, zero at the critical point. Generally, the 
thermodynamic factor displays a well pronounced minimum 
near the Widom line of the studied mixtures, even with only 
0.5 mol% of the solute present. For example, CO2 + 
naphthalene reaches a thermodynamic factor value as low as 
T = 0.45 at T = 320 K along p = 9 MPa. This minimum can 
be related to large density fluctuations and the presence of 
clustering effects near the Widom line. 
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Figure 1: Density dependence of the Fick diffusion coefficient of 
benzene diluted on CO. Molecular simulation data for the Fick 
diffusion coefficient along the isobar p = 9 MPa are compared with 
experimental data and predictive equations (Wilke et al. 1955, 
Funazukuri et al. 2006). 


Fig. 1 shows the density dependence of the sampled Fick 
diffusion coefficient along the p = 9 MPa isobar for the 
supercritical CO2 mixed with benzene in comparison to the 
available experimental data at the same density, but at 
different temperatures and pressures. The anomalous slow- 
down of Fick diffusion near the Widom line observed through 
that experimental data is qualitatively replicated by 
simulation at a finite mole fraction. The observed diffusive 
slow-down becomes less pronounced with decreasing solute 
mole fraction and disappears at the infinite dilution limit, as 
the thermodynamic factor reaches the value of unity. 


Molecular simulation data show a clear breakdown of the 
Stokes-Einstein relation, i.e. decoupling of diffusion and 
viscosity around the Widom temperature, where the transition 
from a high- to low-density fluid occurs, cf. Fig. 2. 
Furthermore, molecular simulation results indicate a strong 
correlation of the Fick diffusion coefficient on the self- 
diffusion coefficient of pure CO2 and the molecular mass of 
the solute, cf. Fig. 3. This relation can be useful for the 
development of an empirical equation to access the Fick 
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diffusion coefficient only based on the properties of pure 
scCOz, the molecular mass of the solute and mixture density. 
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Figure 2: Dependence of the infinite dilution Fick diffusion 
coefficient on the self-diffusion coefficient of pure CO2 for methane, 
ethane, isobutane, benzene, toluene and naphthalene diluted in 
scCOz. Open symbols represent molecular simulation data along the 
9 MPa isobar and temperatures between 290 K and 345 K. The lines 
are linear fits to the data. 
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Figure 3: Ratio of the intra-diffusion coefficient of the solute to the 
solvent for COz2 mixtures with 0.5 mol% of each benzene (red), 
toluene (blue) or naphthalene (green). Solid and open symbols 
represent the simulation data along the isobars p = 9 and 12 MPa, 
respectively. 


Thermodiffusion 


As a second-order effect, the mass flux induced by the Soret 
effect is relatively small and therefore challenging to evaluate 
with equilibrium molecular dynamics simulation. 
Nonetheless, this technique allows for the determination of 
the different contributions to the thermal diffusion coefficient, 
and thus to the understanding of the associated microscopic 
mechanisms. The Soret coefficient relates the thermal 
diffusion coefficient D’ to the Fick diffusion coefficient D,; i 
S™=D"/Dj,;. Therefore, it is expected that the Soret 
coefficient shows an anomalous behavior in the near-critical 
region, as observed for the Fick diffusion coefficient. 
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First molecular simulation results for some of the studied CO2 
mixtures support this assumption with the sampled Soret 
coefficient showing a strong temperature dependence and 
anomalous behavior near the Widom line. Further, analysis of 
the different terms contributing to the thermal diffusion 
coefficient identified the enthalpic term as the main 
contribution. Thus, an accurate assessment of the partial 
molar enthalpies is essential for the calculation of the Soret 
coefficient employing equilibrium methods. 


Conclusions 


Molecular dynamics simulation techniques were employed 
for the determination and analysis of the Fick diffusion 
coefficient of binary mixtures of supercritical CO and 
different solutes. In the near-critical region, it was found, even 
for highly diluted mixtures, that the thermodynamic factor 
may strongly differ from unity. These non-idealities induce 
an anomalous slow-down of the Fick diffusion coefficient in 
this region. Similarly, first simulation results suggest an 
anomalous behavior of the Soret coefficient in the near- 
critical region. 
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Introduction 


Thermo-diffusion is a mass transport phenomenon which 
appears in multi-component mixtures subjected to a thermal 
gradient. This phenomenon is also called the Soret or Ludwig- 
Soret effect (1856). This effect leads to species separation of 
the fluid mixture. It was the subject of numerous works. This 
mass transport takes place towards the hot or cold walls, 
depending on the species involved, their thermodynamic 
characteristics and the thermal boundary conditions. The 
presence of convection and the Soret effect lead to another 
phenomenon called thermo-gravitational diffusion, resulting 
in greater species separation compared to the one obtained in 
pure thermodiffusion. 

To increase the species separation, Khouzam et al. (2013) 
used mixed convection. Sioud et al. (2022) showed that the 
species separation obtained in spherical annular columns is 
more important than in the usual vertical thermogravitational 
column (TGC). Mojtabi et al. (2019) studied a new 
configuration of the horizontal porous layer by moving the 
cavity walls at constant and opposite velocity. Seta et al. 
(2020) studied the species separation stability for different 
mixtures with negative Soret coefficient, in a parallelepipedic 
thermogravitational microcolumn. 

The present study deals with forced convection imposed to a 
binary fluid saturating a porous medium and flowing in a 
cavity whose vertical plates are maintained at constant and 
different temperatures. 

This coupling leads to species separation continuously. This 
process could be used in multiple academic and industrial 
applications. 


1. Mathematical formulation 

In the present study, we consider a water-ethanol mixture 
(60.88% water) flowing through a parallelepipedic cavity, of 
L=2.5 mlength, e = 6-107?m thickness and H = 0,314 
mheight with e << H << L, saturating a porous media of 
a porosity e* = 0.4 Both lateral vertical walls are 
isothermal and maintained different 
temperatures as seen in Fig.1. 


respectively at 


Figure 1: Schematic drawing of the studied cell. 


This mixture is incompressible. We neglect viscous 
dissipation and Dufour effect. 
The Boussinesq approximation is assumed to be valid; in the 
convection generating term, the density varies linearly with 
the local temperature and the local mass fraction. 
The mathematical formulation of the physical problem is 
given by: 

V.V=0 


V= ~ (VP — pol1 — Br (T — To) — Bc (C — Cy) 1g) 
(pc) or + (pc) pV.VT = A* V°T 


Ge +V. vc) = (D* V2C + Dr*Cy(1 — Cy) V2T) 
where D* and D,* are, respectively, the mass diffusion 
coefficient and the thermodiffusion coefficient of the denser 
component in porous media, V is the velocity of the fluid, T 
is the temperature and C is the mass fraction of the 
constituent of interest. We suppose that C(1 —C) = C)(1 — 
Co). 
In these equations, (pc)* and A* are, respectively, the 
effective heat capacity and the thermal conductivity of the 
porous medium, e*its porosity and K its permeability. yw 
and Py, are the dynamic viscosity and the density, at T = Ty 
of the saturating fuid. 
The boundary conditions of the study are written as follows: 
For velocity 

y=0,e;V.n=0 

z=0,H;V.n=0 

x=0;V.x = Ugep 


For temprerature 


x=0;T=T) 
x= h5 oT 
z=0,H; 7-29 
x=0,C=C, 
For mass fraction x=L; <= 


~ OC * aT 
y = 0,e; D at Pr Co(1 a C5 = 0 
2. Analytical solution 
In order to solve the corresponding problem analytically, 


when the flow reaches the massic established regime, the 
forgotten effect hypothesis (Platten et al. 2003) and parallel 
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fow approximation (Mojtabi et al. 2019) were taken into 
consideration. 
V = w(y)z + Uden 
P= fo) 
C=mz+h(y) 
m designates the mass fraction gradient along z. 


Teota — T; 
T= ocr hot y+ Th ot 


wy) — KgBr (cota = Thot)(2Y =: e)/24ev 
(Teaté | Trot) mK gBr(e* = 6ey* + 4y), 


nS D* 24eu 
Dr*C(1—C 1¢ ) 
Tt 0 0 e 2 


_ —10K guBrDr* CoA — Co) (Teta — Tho)? 
(KegPr(T cota -_ Teed? + 120(@vD*)? 


3. Results and discussion 


Numerical simulations are carried out with a finite elements 
method using COMSOL Multiphysics software and the 
results are compared to the analytical ones as seen for vertical 
velocity in fig.2 and mass fraction in fig.3. The evolution of 
the mass fraction at the bottom and top of the column as a 
function of time in x=L is given in fig.4 

We observe a good agreement between numerical and 
analytical results in 2D. 3D numerical simulations are being 
developed to validate used assumptions. 


rf 0.002 0.004 0.006 


Figure 2: Numerical (dots) and analytical (line) of the vertical 
component of the velocity function of y 


Figure 3: Numerical (dots) and analytical (line) mass fraction profile 
as a function of z for x = L 
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Figure 4: Evolution of the mass fraction at the bottom and top of the 
column as a function of time in x=L. 


Conclusions 


This work presents a new procedure allowing the species 
separation in a continuous way. For this, the binary fluid 
initially at a uniform concentration is introduced with a 
constant velocity through the vertical slot of the 
thermogravitational column at x = 0 and comes out through 
the opposite slot at x = L, once mass fraction equilibrium is 
reached. 

It is possible to juxtapose and connect several similar 
columns to achieve complete species separation. 
Also, the fact of using porous thermogravitional columns, 
makes it possible to significantly increase the quantity of 
separated product since the thickness e of the porous cavity 
can be 6 to 7 times greater than the one of a cavity with the 
binary fluid alone. 
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Introduction 


The Enskog solutions to the Boltzmann equations have for 
many years provided a route to predicting the effects of 
thermal diffusion. By using more complex interaction 
potentials, predictions for dilute gases have become 
increasingly accurate for real fluids. However, for dense 
gases, predictions have thus far been limited to the use of the 
hard sphere potential. 


The Revised Enskog Theory (RET) for multicomponent 
mixtures, developed by van Bejren and Ernst, uses the radial 
distribution function (rdf) at contact to account for the effects 
of increased density in hard sphere mixtures. By using 
Barker-Henderson perturbation theory to calculate the rdf at 
contact, we extend the Revised Enskog Theory to Mie fluids 
(RET-Mie). 


Several equations of state have been developed on the basis 
of the Mie potential, and the parameters of the interaction 
potential have been regressed to pure-component equilibrium 
properties. Using these parameters as input, RET-Mie is fully 
predictive for transport properties. 


Comparing the model to an extensive set of gas phase 
experiemental data (Vargaftik, 1975, Oost et al. 1969/1972), 
it is found that predicted thermal diffusion ratios are typically 
within 10-20 % of the reported measurements. It is also found 
that the predicted thermal diffusion factor differs by less than 
15 % from simulation results for Lennard-Jones mixtures, 
even at densities significantly exceeding the critical density. 


Theory 

The Boltzmann equations for a multicomponent system 
consisting of particles with non-negligible covolume may, in 
the absence of external forces, be written as, (de Haro et al. 


1983) 
(= + =f = > JuGofi: 
j 


where u; is the velocity of a particle of species i, f; is the 
velocity distribution function (vdf) of species 7, and 


y= [ff x@+ alone ov) 


— x(r — oy) fi(r — oi; ) fj (r)bdbdedu; 


is the streaming operator, which describes the effect of 
collisions on the vdfs. Here, xy is the radial distribution 
function, oj; is a length associated with the size of the 
particles, for hard spheres taken to be the particle diameter, 


and b and € are the collision coordinates, as illustrated in 
Fig. 1. 


dbde 


Figure 1: The coordinates describing a binary collision. 


An approximate solution to this set of equations is found by 
means of the Enskog solution method. From these solutions, 

it becomes apparent that the “contact diameter”, o;;, plays 
two distinct roles: That of modifying the frequency of 
collisions, through the rdf at contact, and that of enhancing 
transport as energy and momentum are transferred across a 
non-zero distance from one particle to the next during 
collisions. 


We employ Barker-Henderson perturbation theory to 
evaluate the rdf at contact for particles interacting through a 
Mie (generalized Lennard-Jones) potential. Regarding the 
instantaneous transport between particles at the moment of 
collision, we propose a manner in which to find a reasonable 
distance that may be considered the “contact diameter” for 
colliding Mie particles. Thus, we may employ the RET to 
evaluate the transport coefficients of a dense, multicomponent 
mixture of Mie particles. 


Mie potential parameters 


All predictions in the current work have been made using 
potential parameters regressed to equilibrium properties. 
Calculation of transport properties thus requires no a priori 
knowledge about the transport properties of the investigated 
systems. 
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Implementation 


The proposed model, RET-Mie, has been implemented as an 
Open-Source package, and may be found on GitHub under 
the ThermoTools project. 


Figure 2: The thermal diffusion factor of several isotopic Lennard- 
Jones mixtures at T* = 2.0. Predictions (lines) using RET-Mie. 
Marks are simulation results (Hoang and Galliero, 2022). 


Results 


The RET-Mie has been tested in its ability to reproduce the 
transport properties of a wide variety of systems. Not only 
thermal diffusion has beem examined, but also diffusion 
coefficients, shear viscosities and thermal conductivities. 


As shown in Fig. 2, the predicted thermal diffusion factor of 
isotopic LJ mixtures is in good agreement with simulation 
results (Hoang and Galliero, 2022), even up to a reduced 
density of p* = 0.7. 


—— He/Ne Ne/Ar 
—t— He/Ar —— Ne/Kr 
—— He/Kr —@ Ar/Kr 


200 400 600 800 


Figure 3: The thermal diffusion factor of several real mixtures, as 
predicted using RET-Mie (lines) and computed using ab initio 
potentials (marks) (Song et al. 2011). 
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Thermal diffusion factors were predicted for several real 
mixtures at temperatures ranging from 100 K to 1000 K, and 
compared to computations using ab initio quantum chemical 
potentials (Song et al. 2011). As shown in Fig. 3, the 
agreement with their computations is excellent. 


Conclusions 


The RET-Mie has so far shown itself to be a reliable model 
for predicting thermal diffusion coefficients, and other 
transport coefficients, without requiring any a_ priori 
knowledge about the transport properties of the system in 
question. 


Several possible extensions, both involving molecules with 
internal energy, molecules with significantly non-spherical 
interaction potentials, and systems with significantly large 
thermal de Broglie wavelengths deserve further attention. 
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Motivation 

Our study of ionic solutes is motivated by the most important 
practical application of thermodiffusion, where it is used to 
monitor protein-ligand reactions. Proteins are complex 
molecules that contain ionic as well as non-ionic groups. 
While non-ionic solutes in water have been extensively 
studied recently (Niether and Wiegand, 2019), ionic solutes' 
concentration and temperature dependence have not been 
investigated systematically. For non-ionic compounds, a 
strong correlation between thermodiffusion and hydration 
was found (Niether and Wiegand, 2019). 


Figure 1: Schematic comparison of the temperature dependence 
of Sr for non-ionic and ionic solutes at different concentrations: 
low (dotted line), intermediate (dashed line), and high (solid line). 


Comparison of non-ionic and ionic solutes 

We found one striking difference between non-ionic and ionic 
solutes looking at the effect of concentration on the 
temperature dependence of the Soret coefficient S;, as 
illustrated in Fig. 1 (Mohanakumar et al. 2021). For a typical 
non-ionic solute in water, the behavior of Sy changes from 
increasing with temperature to decreasing with temperature 
as the concentration increases. This is correlated with the 
hydration of the solutes which decreases as the concentration 
increases. Only very hydrophilic non-ionic solutes have Sr 
values that increase with temperature for all concentrations. 
In contrast, the Soret coefficients of ionic solutes show the 
typical temperature dependence of very hydrophilic solutes 
over the entire concentration range. For salts with a high 
degree of dissociation we might have a tightly bound first 
hydration layer, which leads to a highly hydrophilic entity. 
For less dissociated salts it might be explained by cluster 
formation of the salts with increasing concentrations. Even at 
high salt concentrations the clusters as a whole are hydrated 
at their surfaces by water, but the total exposure to water is 


less as the surface decreases when more ions are part of larger 
clusters. 
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Figure 2: (a) S'r values of all studied systems plotted as function 


of log P. Note, that log P is the sum of an ionic and non-ionic 
contribution. First order polynomials of concentration and 
temperature have been used to fit the data using Eq.1. (b) 
Sequence of the anions based on S'r for the two investigated 
cations in comparison with the Hofmeister series. 


Anion and Cation influence on Sr 

We investigated systematically the concentration and 
temperature dependence of the thermodiffusion of aqueous 
solutions of various potassium and_ sodium _ salts 
(Mohanakumar et al. 2021, Mohanakumar et al. 2022a). 

To describe the temperature and concentration dependence 
we used an empirical Ansatz suggested by Wittko and Kohler 
(Wittko and Kohler, 2007) 


Sp(m,T) = a(m)B(T) + Sp (1) 
With polynomial serial expansions for a(m) and B(T) 


a(m) = aym+a,m? +a,m3 +-: 
Bm) = 1+5),(T —T.) + bo(T — Tp)? + (2) 
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m is the molality, Ty is an arbitrary reference temperature, 
set to To = 25°C and St is a temperature and concentration 
independent constant. Note, that we set aj =0 as it is 
strongly coupled to St. 

We could describe the temperature and concentration 
dependence of S; of various potassium and sodium salts in 
water using Eq.(1). In Figure 2(a) we display the adjustable 
parameter $4 as function oflog P, with P being the ratio 
of the equilibrium concentration of the solute (salt) in octanol 
and in water. So, a negative logP signifies stronger 
hydrophilicity. We find for all investigated salts a linear 
correlation between Si and log P. This implies that also, for 
ionic solutes, hydrophilicity plays an important role. If we 
compare with the hydrophilicity scale of the Hofmeister 
series, we find a good agreement except for the thiocyanate 
anion, which should be, according to Hofmeister, the most 
hydrophobic anion. In brief, we can state that the hydration of 
the ions plays a significant role. 


Overlapping hydration shells 

Various salts in water exhibit non-monotonic variations of the 
Soret coefficient S; with concentration, which is not 
understood on a microscopic level. We investigated the 
thermodiffusive properties of aqueous solutions of sodium 
iodide, potasstum iodide and lithium iodide, using thermal 
diffusion forced Rayleigh scattering in a concentration range 
of 0.5 — 4 mole per kg of solvent and a temperature range of 
15 to 45°C (Mohanakumar et al. 2022b). In all three cases Sy 
has a minimum at Myj, = 1 mole per kg of solvent. We 
develop an intuitive picture in which the relevant objects are 
the fully hydrated salt molecules (FHP), including all water 
molecules that behave differently from bulk water. Our 
hypothesis is that these FHPs form a random close packing at 
Mmin, Which implies that the outer hydration shell start to 
touch as indicated in Figure 3. Preliminary, somewhat sketchy 
calculations indicate that indeed Soret coefficients begin to 
rise beyond mypjn. Indications are given as to why the model 
will fail at large concentrations. 


: oe 


molality Min 


Figure 3: Hydrated salt molecules overlapping with increasing 
concentration. The green-red sphere represents the bare salt 
molecule, after adding the blue shell of strongly attached water 
molecules we get the salt particle (HSP), while after adding next 
the outer light blue shell of perturbed water we arrive at the 
hydrated salt molecule, called FHP. At concentrations above mmin 
the outer shells overlap as shown on the right side. 
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Conclusions 

We have studied the thermophoretic properties of various 
salts in water over a range of temperatures and concentrations. 
Although the temperature dependence of the Soret coefficient 
of the ionic compounds does not change in the same 
pronounced way as function of their hydrophilicity observed 
for non-ionic solutes, we find a linear correlation between Si 
and log P. Most likely, the hydration shell of ionic solutes is 
more tightly bound to the ions than in the case of non-ionic 
solutes, so we find a similar temperature dependence of the 
Soret coefficient for all concentrations. Additionally, 
overlapping of the hydration shells might also be responsible 
for the occurrence of a minimum of Sy with concentration. 
However, this hypothesis needs to be quantified by computer 
simulations. 


Acknowledgements 


We thank Fernando Bresme, Jan Dhont and Jutta Luettmer- 
Strathmann for fruitful and helpful discussions. 


References 


D. Niether, S. Wiegand, Thermophoresis of biological and 
biocompatible compounds in aqueous solution, J. Phys. 
Condens. Matter, 31, 503003 (2019). 


S. Mohanakumar, J. Luettmer-Strathmann, S. Wiegand, 
Thermodiffusion of aqueous solutions of various potassium 
salts, J. Chem. Phys., 154, 84506 (2021). 

S. Mohanakumar, S. Wiegand: Towards understanding 


specific ion effects in aqueous media using thermodiffusion 
The Eur. Phys. J. E 45(2), 10 (2022a). 

G. Wittko, W. Kohler, On the temperature dependence of 
thermal diffusion of liquid mixtures, Europhys. Lett. 78, 
46007 (2007). 

S. Mohanakumar, H. Kriegs, W. J. Briels, S. Wiegand, 
Overlapping hydration shells in salt solutions causing non- 


monotonic Soret coefficients with varying concentration 
PCCP 24, 27380 (2022b). 


36 


15" International Meeting on Thermodiffusion (IMT15) 


Tarragona (Spain), May 2023 


Modelling thermodiffusion in aqueous sodium chloride solutions — best performing 
water models for predicting inversion temperatures 


Alice J. Hutchinson!”, Juan F. Torres!, Ben Corry” 


School of Engineering, Australian National University, Canberra, Australia, felipe.torres@anu.edu.au, 
Research School of Biology, Australian National University, Canberra, Australia, ben.corry@anu.edu.au 


Introduction 

Thermodiffusion is the migration of a species due to a 
temperature gradient and is the driving phenomenon in many 
applications ranging from early cancer detection to uranium 
enrichment. Molecular dynamics (MD) simulations can be a 
useful tool for exploring the rather complex thermodiffusive 
behaviour of species, such as proteins and _ ions. 
Thermodiffusive behaviour of a species can be described by 
the Soret coefficient St, which indicates the magnitude and 
direction of species migration under a temperature gradient, 
and the inversion temperature 7o, which is the temperature at 
which a species switches from displaying thermophilic to 
thermophobic behaviour. Current MD _ models of 
thermodiffusion in aqueous ionic solutions struggle to 
quantitatively predict St and 7o. In this work, we aim to 
improve the accuracy of MD thermodiffusion models by 
assessing how well different water models can recreate 
thermodiffusion in a benchmark aqueous NaCl solution. 
Beyond this, we aim to use the best performing water model 
to explore how the inversion temperature 7o of a NaCl 
solution changes with concentration of the solution. 


Which water model is best? 

To improve predictions of thermodiffusion, the accuracy 
limitations of MD modelling must be considered. MD is 
based on classical mechanics and, therefore, inherently 
approximates real-world molecules and their behaviour. In an 
aqueous’ electrolyte MD _— simulation, the major 
approximations lie upon the choice of the ion and water 
models (force fields), which define the structure and 
interactions of these molecules. While the exact mechanism 
of thermodiffusion is not well understood, multiple sources 
emphasise the importance of hydrogen bond structures, 
solute—water interactions, and hydration entropy in defining 
thermodiffusive behaviour (Duhr & Braun 2006, Di Leece et 
al. 2017, Niether & Wiegand 2019). The water model and its 
associated ion parameters define the molecule structure and 
the water—water, ion—ion, and water—ion interactions in the 
solution and, hence, are integral to accurately recreating 
thermodiffusion. 


In this work, the thermodiffusive behaviour of three- and 
four-point models from two of the best performing water 
model families were tested (Figure 1). TIPnP-FB 
approximates water from on experimental reference data 
(Wang et al. 2014) and OPCn approximates water from 
quantum mechanically calculated reference data (Izadi et al. 
2014, Izadi et al. 2016). We also tested two commonly used 
models: TIP3P (Huang et al. 2013), an old but regularly used 
water model that is expected to perform more poorly, and 
SPC/E (Berendsen et al. 1987), the three-point model 
commonly used in other MD thermodiffusion works (R6mer 
et al. 2013, Di Leece et al. 2017). Specifically, we sought to 
compare the performance of the chosen water models in 


reproducing 70 of 0.5M aqueous NaCl solutions and St of 0.5, 
2, and 4M aqueous NaCl solutions, as compared to the 
experimental values obtained by Romer et al. (2013). 


We directly observe the thermally driven diffusion of the ions 
in solution, by developing a quasi-linear temperature profile 
over the simulation. The temperature profile is achieved via 
two spatially defined thermostats controlled to different 
temperatures (one hot, one cold). Molecules diffuse across the 
temperature profile to reach a steady state distribution where 
thermally driven diffusion is opposed by Fickian diffusion. To 
are found by a maximum in the steady state concentration 
profile of the ions. The Soret coefficient is derived by fitting 
the concentration profile versus temperature. Each water 
model predicted a different ion distribution yielding different 
To and Sr (Figure 2). By comparing the modelled St with 
published experimental values, we determined TIP3P-FB to 
be the water model that best recreates thermodiffusion in 
aqueous NaCl solutions (Hutchinson et al. 2022). 


ie, {_.. SPCE 


@ Sodium ion © Chloride ion 


Temperature gradient 


Initial ion distribution 


sy Ree 


Steady-state ion distribution Jt 
= = 
i P.O 


Localtemperature: T<T To T>T 
lon behaviour: Thermophilic ————» <———Thermophobic 
Soret coefficient: S7< 0 $7 =0 S7>0 


Figure 1: Zop: Depiction of MD model of aqueous NaCl solution. 
The model volume is occupied by water molecules (light grey), 
sodium ions (yellow), and chloride ions (purple). Six different tested 
water models are represented on the right hand side by different 
colours. A cold and hot thermostat (highlighted blue and red areas, 
respectively) are spaced symmetrically within the simulation to 
create a linear temperature profile over the system. Ions diffuse 
within the temperature profile. Bottom: Depiction of thermophilic 
and thermophobic behaviour of an aqueous NaCl solution under a 
linear temperature profile that encompasses an _ inversion 
temperature. The direction of thermally driven ion diffusion is 
indicated by the sign of the Soret coefficient, S;, and depends on 
whether the local temperature of the solution, T, is greater or less 
than the inversion temperature, To. 
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Figure 2: MD results for Soret coefficient and inversion temperature 
in 0.5M NaCl solutions as a function of temperature for each tested 
water model (Hutchinson et al. 2022). Experimental data from 
Romer et al. (2014) is shown in black for the comparable 
concentration of 0.5 mol kg! NaCl. A measure of distance between 
the predictions of each water model and the experimental curve is 
quantified by a y? value in the right-hand side table, where smaller 
values imply the water model predictions are closer to the curve. 


How does the inversion temperature of an aqueous NaCl 
solution depend on ion concentration? 
Following the results in Hutchinson et al. (2022) outlined 
above, the same general MD simulation set-up is used to 
explore how 7» of aqueous NaCl solutions depend on ion 
concentration in the range of 0.05—2.0M. 


It has been proposed that the Soret coefficient is determined 
by two competing, temperature dependent values — the 
entropy of hydration and the entropy of ionic shielding (Duhr 
and Braun, 2006). The entropy of hydration is defined for a 
single ion in solution. Our preliminary MD simulations 
suggest that hydration entropy alone cannot cause inversion. 
A competing entropy of ionic shielding is required for 7o to 
be defined, a value that is not defined for single ions in 
solution. Therefore, we hypothesise that in the infinite 
dilution limit (concentration,c > 0), solutions can only 
behave thermophobically (i.e. Zo is not defined). 


Initial MD modelling of NaCl aqueous solutions of 
concentrations 0.1, 0.2, 0.5, 1, and 2M supported this 
hypothesis with a seemingly asymptotic drop in 7o at the 
lowest concentration of 0.1M (Figure 3). The preliminary 
results reproduce the known 7o at moderate concentrations, 
and the known reduction in inversion temperature with 
increasing concentration of NaCl. However, the simulation 
error was too large to be certain of these results and hence 
further simulations are to be undertaken. 


Conclusions 

All tested rigid non-polarisable water models qualitatively 
reproduced thermodiffusive behaviour. Simulations using 
TIP3P-FB, followed by TIP4P-FB and OPC3, resulted in the 
most accurate quantitative predictions of thermodiffusion for 
a 0.5M NaCl solution. More accurate thermodiffusion 
properties were obtained by using combined water and ion 
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Figure 3: Inversion temperature of aqueous NaCl of varying ion 
concentrations predicted by MD simulations. Measured inversion 
temperatures (black dots) are displayed as an average with standard 
deviation amongst replicates of the simulation. Red line is drawn as 
a guide to the eye, supporting the hypothesis of a lack of inversion 
in single solute solutions and proposing a defined thermophobic and 
thermophilic behavioural region for NaCl solutions. 


models that accurately recreate standard properties of the ion— 
water interactions, such as hydration free energy and ion— 
oxygen distance, rather than standard properties of water. 


The results provide further confidence that MD simulations 
can capture molecular scale details of the interactions 
between ions and water molecules to achieve a better 
understanding and predictions of thermodiffusion. We further 
applied such MD simulations to explore the relationship 
between inversion temperature and ion concentration in NaCl 
aqueous solutions. Preliminary results support the hypothesis 
that single solute solutions do not exhibit inversion behaviour 
and that complex thermodiffusion of NaCl solutions, and 
likely other aqueous solutions, cannot be described by single 
solute properties alone. 
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Introduction 


Einstein-Helfand (EH) and Green-Kubo (GK) formulas [1,2] 
have been widely used to calculate the values of macroscopic 
transport coefficients from equilibrium molecular dynamics 
simulations [3]. Transport coefficients, such as diffusion 
coefficients, shear viscosities or thermal conductivities, are 
obtained from expressions containing time correlation 
functions of microscopic fluxes (e.g. particle flux, stress 
tensor or heat flow), which are sampled from the dynamic 
simulations. Compared to non-equilibrium simulations, 
which reproduce the experimental setup in which e.g. a 
temperature gradient produces a heat flux across the sample, 
EH and GK are calculated from equilibrium simulations. The 
EH and GK approaches do not require specific boundary 
conditions to create the density, velocity or temperature 
gradients suitable to generate the conjugated flow. Therefore, 
one single equilibrium simulation can provide a set of 
transport coefficients. In this presentation we show that the 
standard GK and EH formulas can be applied to Dissipative 
Particle Dynamics (DPD) models, whose dynamics is 
expressed in terms of Langevin-like equations, containing 
dissipative and random contributions, contrary to the general 
belief. 


Shear viscosity 


Despite its generalised use in standard molecular dynamics 
simulations, concerns have been expressed about whether the 
presence of irreversible and random interactions forbids the 
use of EH and GK formulas in DPD. For such systems with 
stochastic and dissipative interactions, Ernst and Brito [4] 
proposed modifications of the GK formulas in which the 
transport coefficient is expressed as a summation of different 
contributions, which include conservative and dissipative 
interactions, but not the random contributions. For the shear 
viscosity, they propose 


fy at (1S) +12.) (1.0) - 
(1) 


1 
1 = No + kpTV 
11,2, (0))) 


The first term is the effect of the random contribution and the 
second term (inside the integral) is affected by a negative sign, 
which is at least striking. Using this approach, Jung and 
Schmid [5] applied the GK relation of Eq. (1) for the 
calculation of viscosity for Dissipative Particle Dynamics 
(DPD) models. 


However, in this presentation, we show that the standard EH 
and GK formulas are in fact valid for DPD-like models which 
satisfy only two fundamental conditions, namely, i) the model 


USA 


preserves the conservation of the given property (in the case 
of the viscosity, DPD models preserve momentum 
conservation) and ii) the random contributions satisfy detailed 
balance. Extensive simulations show that the standard 
expressions are consistent if the appropriate form of the fluxes 
is introduced in the formulas, i.e. the dissipative and random 
contributions are also included. For the sake of comparison 
with Eq. (1), we can prove that the expression 


= a fy dt (Il,,(t)M,2(0)) = No + 
1 


sive tH (MEO + 12©)(EO) + 12,0) + 
11,8(0))) (2) 


is theoretically equivalent to Eq. (1). In Fig. (1) we show the 
numerical comparison between the two expressions. 
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Figure 1: Comparison between values of the viscosity 
obtained from our proposed EH formula and non-equilibrium 
simulation results of Figure 5 in Ref. [5]. Blue circles 
correspond to simulations performed with the isothermal 
DPD with y = 5, T = 1, and the same time step as in the 
original reference, ot = 0.01. The errors are of the size of the 
symbols. Green triangles stand for the values obtained from 
GenDPDE simulations with large heat capacity Cv = 1000 
and an interparticle thermal conductivity k = 50. 


Thermal conductivity 


DPD with energy conservation (DPDE) [6-9] is a particular 
type of DPD model that includes an internal energy and 
particle temperature along with conservative, dissipative and 
random interactions. DPDE, and its generalisation to complex 
systems (GenDPDE [10-12]), is of particular interest in 
mesoscopic fluid simulations and has applications in several 
areas of interest including fluids with chemical reactions, 
multiphase flow and heat transport under complex conditions. 
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Following the same line of reasoning, we present the 
derivation of the formulas for the EH and GK relations for the 
thermal conductivity (Figure 2) in DPDE and GenDPDE. The 
calculation of the transport coefficients from equilibrium 
simulations was also contrasted with non-equilibrium 
calculations [13]. We observe very good agreement of the 
proposed equilibrium EH and GK relations under different 
parameter ranges, proving the validity of the expressions. In 
overall, we present an extensive analysis of the validity of EH 
and GK formulas in DPD models and we extend its 
application to DPDE and GenDPDE models. 
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Figure 2: Thermal conductivity as a function of the density 
for values of the particle thermal conductivity k = 5 and 50. 
The circles represent the results obtained from non- 
equilibrium DPDE simulations (blue circles « = 5 and red 
circles k = 50), the triangles from equilibrium DPDE 
simulations using our proposed FH relation (grey triangles « 
= 5 and green triangles « = 50), and the predicted theoretical 
value from mean-field approximation (dotted black line). 


Conclusions 


Our theoretical analysis, together with the broad range of 
parameters studied, are a strong case against the false belief 
that DPD-like systems cannot be treated on the same grounds 
as conservative systems. The random contributions in the 
dynamics introduce some peculiar behaviours, which are not 
present in conservative systems. However, the general 
principles of conservation and detailed balance suffice to set 
the appropriate framework for the analysis of transport in 
dissipative and stochastic systems, with no ambiguity. 
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Introduction 

Thermodiffusion may induce non-homogenous concentration 
distributions due to temperature gradients. Despite being 
discovered more than 150 years ago and having a rich history 
of observation, no universal explanation has been found for 
the physics behind thermodiffusion. For many species, 
including ions (R6émer et al., 2013), colloids (Putnam, Cahill 
and Wong, 2007), polymers (Qian et al., 2018), and 
biomolecules (Duhr and Braun, 2006), the Soret coefficient 
St can either be positive (thermophobic) or negative 
(thermophilic). The temperature at which species changes 
thermophobicity (i.e. migration direction) is termed the 
inversion temperature 7o and its value depends on the 
concentration and temperature. Note that To is generally low 
for most species in aqueous solutions, e.g. < 25°C in salt 
solutions and polymers (Putnam, Cahill and Wong, 2007; 
Romer ef al., 2013). To occurs due to the interaction between 
the species and effects of its thermal environment (De Miguel 
and Rubi, 2019). Data on Zo could help understand the 
underlying physics of thermodiffusion. However, measuring 
To may be challenging because the required low temperature 
control, negligible thermodiffusive transport around 7o, i.e. 
the measurable AC is very small (Romer ef al., 2013), and 
limitations in some optical methods such as beam deflection 
which only give limited one-dimensional information 
(Kolodner, Williams and Moe, 1988). Here, we demonstrate 
a fast and sensitive digital interferometric technique that can 
measure 7o in a transparent aqueous solution. Furthermore, St 
can be extracted for a temperature range at a given 
concentration from a single experiment using this technique. 


Proposed experimental methods 

A transparent binary or multi-component solution is placed 
inside a Soret cell (Figure 1) whose upper and lower 
boundary temperatures are controlled with Peltier modules, 
thermistors and a proportional-integral-derivative control 
system. A temperature stability of 3 mK is achieved. The heat 
extracted by the Peltier modules when cooling the Soret cell 
is effectively removed with a water circulation system. 


Phase-shifting interferometry (PSI) is a temporal digital 
interferometry technique that can resolve the phase difference 


Water Nickel- 
cooling plated Sealant 
blocks copper 


Quartz Peltier 
cell modules 


Figure 2: Soret cell. Green dashed arrows indicate the 
embedded temperature sensor in the copper blocks. 


Ay between the test and the reference beams; see more details 
in our recent work (Torres et al., 2020). In the experiments, 
phase-shifted data (as those shown in Figure 2) is obtained 
from three interferograms acquired within 100 ms, which is 
much shorter than the time taken for the temperature and 
concentration fields to fully develop in a 4 mm high Soret cell. 
The experimental procedure is shown in Figure 2. First, we 
set a no-fringe condition for an isothermal and _ iso- 
concentration solution, as shown in (a). Since the Lewis 
number is in the order of 100, the fringes that appeared in the 
first 1 min since applying a temperature difference, e.g. AT of 
5 K in (b), can be suppressed by adjusting a test beam mirror 
(Torres et al., 2013), as shown in (c) for a quasi-homogeneous 
solution having a linear temperature profile. After 30 min, the 
concentration field has fully developed, such that phase 
differences in (d) are due to concentration gradients only. 


In the unwrapped phase difference (i.e. phase map), the phase 
y for all pixels is available and the horizontal average yields 
the phase map as a function of height, w(v). The phase map 
just after adjusting the mirror is yo (Figure 2c) and after the 
concentration field has fully developed is wena (d). 
Furthermore, in a convectionless environment with a linear 
temperature gradient across the cell height, we consider that 
there is only mass flux in the vertical direction, meaning that 
the horizontal mass flux is neglected within the field of view 
that captures the top and bottom boundaries. Note that y is 
directly related to the vertical concentration profile, which 
satisfies conservation of mass, i.e. fr Wo = ti Wena. The 
mean vertical values of yo and wena are then subtracted from 


Adjust mirror | 


(c) 


Figure 1: Experimental procedures. (a) First, Tiop and Tyim are set 
to Tmean thus homogenous temperature and concentration is 
achieved in the solution. The test beam is adjusted to set a no-fringe 
condition. (b) Tiop and Tytm are set to TmeantAT and Tmean—AT, 
respectively. Fringes start to appear quickly within the Soret cell 
due to temperature gradients. The development is less than 1 min. 
(c) Once the temperature field has reached steady state, the test 
beam is quickly adjusted to remove all fringes. (d) After 30 min 
when the cell height H is 4 mm, the thermodiffusive concentration 
field becomes fully develop. Fringes/gradients in the unwrapped 
phase-shifted data are the result of the concentration field. 
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-e 3: yas a function of y and 7 for NaCl concentration 0.5 M. 
ae mean values of yo and wena are subtracted from the yo and 
espectively to obtain yo and w’ena. (b) Ay’ = wend — wo. Grey 
‘s indicate the direction of NaCl migration and 0Ay’/ch is 
proportional to the magnitude of mass flux. Between 12 °C (Trop) 
and 7 °C (Totm) NaCl shows thermophobic behaviour. (c) When Trop 
= 15°C and Titm = 10°° 
thermophilic behaviour. Direct inspection yields 7 ~ 11.3 °C. 


yo(v) and wend(y), respectively, to reflect mass conservation 
(y’). This process is depicted in Figure 3a. Since the contrast 
factor can be assumed to be constant in a small concentration 
range (expected in a thermodiffusion experiment with St ~ 0), 
the obtained Aw’ is a direct measurement of the local 
concentration change. In addition, with the assumption of 
linear temperature gradient across the channel height, the 
vertical coordinate y can be transformed into temperature T 
by interpolation using the top Tiop and bottom btm 
temperatures. In the experiment of Figure 3b, based on 
literature values (Caldwell, 1973; Rémer et al., 2013), we first 
estimated that 7o for 0.5 M NaCl aqueous solution should be 
between 7 and 12 °C and used these temperature as 7otm and 
Twop, respectively. However, NaCl demonstrated thermophilic 
behaviour when 7’< 11 °C as 0Aw’/oh is positive (b). Between 
11 and 12 °C, dAw’/oh is close to zero, suggesting near-zero 
St, i.e. To is around this location. In (ce), the boundary 
temperatures are re-set to confirm the previous observation. A 
maximum Ay’ is obtained at the location corresponding to T 
= 11.3 °C. Above 11.3 °C, NaCl is thermophobic (OAy’/oh > 
0), and below 11.3 °C NaCl is thermophilic (OAy’/oh < 0). 
With known contrast factors, C() can be obtained and based 
on a steady-state condition, the local Sr (i.e. at a given 


temperature and location) can be calculated as a 
C(1-C)dT 


Thus, the dependence of St on 7' can be precisely measured. 


Results 

Following the procedures outlined in Figure 2 and data 
processing steps in Figure 3, 7o as a function of concentration 
was measured for both aqueous NaCl and KCl from multiple 
experiments, with the results shown in Figure 4. We can 
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for NaCl and 0.5 M for KCI. Beyond these values, To starts 
dropping slowly with increasing concentration. Experimental 
data on 7o for the investigated solutions is only available for 
NaCl/H20 at 0.5 M in literature (Caldwell, 1973) for which a 
good agreement was obtained. S1(7) was also extracted from 
~ concentration profile obtained using PSI. 


aclusions 

‘ough PSI and post-data processing, a steady concentration 

file in a Soret cell can be extracted with high accuracy. 

asured Jo for both investigated binary solutions seems to 
follow an interesting common trend with increasing C: 
beginning with a sharp increase, reaching a maximum 7o and 
then followed by a gradual decrease. Further work on 
quantifying errors in this measurement is required. 
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Introduction 


The thermodiffusion of ions in liquids is difficult to predict 
due to an an incomplete understanding of its microscopic 
origins. It is widely hypothesised that thermodiffusion is 
driven by changes in solvation and ionic shielding entropy, 
creating an entropic force (Duhr & Braun, 2006; Niether et 
al., 2017). Duhr and Braun (2006) suggest that large negative 
ionic shielding entropies account for thermophobicity at high 
temperatures, whilst large positive hydration entropies 
explain thermophilicity at low temperatures. In the case of 
aqueous solutions, structural changes in the hydrogen bond 
network are also thought to contribute to the direction of 
thermodiffusion (Niether et al., 2017). It is argued that 
structure formation at low temperatures minimises the free 
energy, whilst structure breaking due to larger translational 
entropic gains minimises the free energy at high temperatures 
(Wang, Kriegs, & Wiegand, 2012). This free energy change 
determines the probability of particle diffusion under 
equilibrium thermodynamics. 


Alternatively, some authors attribute the mechanism of 
thermodiffusion of molecules in ionic solution to the creation 
of local and global electric fields, reducing this to an 
electrophoretic effect (A. Putnam & G. Cahill, 2005; Reichl, 
Herzog, Gotz, & Braun, 2014; Wiirger, 2008). In particular, 
Reichl et al. (2014) propose extra contributions to 
thermodiffusion arising from the Seebeck effect and the 
capacitor model of Debye screening. Crucially, their model 
suggests thermophoresis is a property of bulk solutions, rather 
than being singularly dependent on molecule-solvent 
interactions. 


Whilst the Soret effect has been extensively studied in 
mixtures and alkali halide solutions, the literature 
conspicuously lacks thermodiffusion experiments involving 
single ionic particles. Such experiments could resolve 
whether thermophoresis only occurs in bulk solutions. 
Further, single-ion studies would remove ionic shielding 
effects, allowing direct studies of the solvation entropy and 
free energy on inversion alone. 


Molecular dynamics (MD) simulations can accurately model 
the atom scale interactions essential to thermophoresis and 
have been shown to reproduce thermophoretic behavior 
(Artola & Rousseau, 2013; Belkin, Chao, Giannetti, & 
Aksimentiev, 2014; Di Lecce, Albrecht, & Bresme, 2017; 
Hutchinson, Torres, & Corry, 2022). Such simulations also 
provide a _ viable avenue to study _ single-particle 
thermophoresis, although this is yet to be pursued. 


This study uses alchemical free energy perturbation (FEP) 
methods to isolate the ion hydration free energy in silico, for 
individual cations and anions. This will indicate if 
thermodiffusive behavior is solely a property of bulk solution. 


Isolation of a single ion in a water box removes ion-ion 
interactions, thereby allowing the study of thermodiffusion 
caused purely by solute-solvent effects. This investigation 
tests whether the entropic force mechanism of 
thermodiffusion is validated for single ionic species over a 
range of temperatures. Significant emphasis is also placed on 
determining whether a thermophilic to thermophobic 
inversion point exists for single particles, a transition 
observed in bulk alkali halide solution (Di Lecce et al., 2017; 
Hutchinson et al., 2022). Investigation of this effect for a 
variety of cations and anions also allows an investigation into 
the fundamental ion properties governing these effects. 


Section 1 — Free Energy Perturbation 


Simulations investigating the free energy of ion hydration 
with respect to temperature and ion identity have been 
performed. The hydration free energy at each respective 
temperature was calculated via alchemical FEP methods. The 
TIP3P-FB forcefield has been used here, as it has been shown 
to be reliable for modelling thermodiffusive behaviour 
(Hutchinson et al., 2022). Figure | shows a comparison of the 
hydration free energy of each ion at each temperature is 
compared to that ion’s hydration free energy at 270K. In all 
case, the ions exhibit thermophobic character, having more 
favourable hydration free energies at lower temperatures in 
this infinitely dilute reigime. By comparison, 0.5M NaCl has 
an inversion temperature at 10°C, becoming thermophilic 
below this temperature (R6mer, Wang, Wiegand, & Bresme, 
2013), something that is reproduced using MD simulations 
with the same parameters (Hutchinson et al., 2022). 
Furthermore, this model system shows that thermodiffusive 
force, given by the slope of the hydration free energy, is ion 
specific (Figure 1), indicating thermodiffusion might be 
utilised in the separation of different ion species. 
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Figure 1: The temperature dependence of the hydration free energy 
of ions. 
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Section 2 — Ion parameters 


The ion-specific nature of this temperature dependence effect 
in the hydration free energy has strong correlations with the 
ion’s radial charge density (Gregory, Wanless, Webber, 
Craig, & Page, 2021) for the monatomic ions presented here 
(Figure 2). This suggests that the thermophobic character 
observed might be predicted by knowing the properties of the 
solute involved, and that the charge density is fundamentally 
important, at least in the case of monatomic ions. 
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Figure 2: The temperature dependence of the free energy of ion 
solvation correlates with the ion’s radial charge density, p. 


Conclusions 


Lone ions at infinite dilution (achievable via simulation but 
difficult to achieve in experiments) appear to exhibit 
thermophobic character at all temperatures investigated here. 
This contrasts to concentrated electrolytes in bulk solution 
that exhibit a temperature inversion, having both thermo- 
phobic and -philic character. Furthermore, the nature of the 
thermophobicity appears to be dependent on the identity of 
the ion (1.e., ion-specific), correlating with the ion’s radial 
charge density. By using simulations, we have been able to 
isolate some of the key factors governing thermophobicity of 
ions in water, providing a framework for predicting the 
thermodiffusion of aqueous ions in dilute aqueous 
environments. 
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Introduction 


The coupling between heat and mass transport at the 
nanoscale is a generic problem found in_ several 
physicochemical scenarios, ranging from molecular motors to 
heterogeneous catalysis, among others, and includes relevant 
chemical engineering applications. The analysis of these 
scenarios may be challenging using molecular dynamics 
simulations, whose computational cost is prohibitive for 
systems with large characteristic time and length scales. 
Coarse-grain (CG) modelling represents a valid alternative 
for such cases where molecular simulations are 
computationally limiting. 


GenDPDE and GenDPDE-M 


Dissipative Particle Dynamics (DPD) was originally 
introduced as a CG Lagrangian method to describe isothermal 
fluid dynamics with thermal fluctuations (Hoogerbrugge et al. 
1992). Heat transport was soon included through energy 
conservation in what is referred to as the energy-conserving 
DPD (DPDE) model (Bonet Avalos et al. 1997). In recent 
years, both methods have been increasingly used in different 
areas, ranging from biology to engineering applications. 

The Generalised Energy-Conserving Dissipative Particle 
Dynamics method, GenDPDE, (Bonet Avalos et al. 2019) has 
been introduced as an extension of DPDE, based on the 
definition of a mesoscopic thermodynamic description at the 
particle level. Within this framework, the particles carry 
internal energy and have variable volumes, which permit the 
definition of a fluctuating, locally-defined temperature and 
pressure. Differences in particle temperatures induce heat 
fluxes between the particles, following Onsager’s Linear 
Theory at the mesoscale. The particle pressure is directly 
related to the forces exerted between the particles, leading to 
temperature-dependent particle potentials that are relevant 
when capturing non-equilibrium phenomena present in 
chemically reacting systems (Lisal et al. 2022) and systems 
undergoing shock compression (Lee et al. 2023). 

Recently, GenDPDE has been further extended to include 
particles of variable composition. The new method, 
GenDPDE-M, (Bonet Avalos et al. 2022; Lisal et al. 2022) 
can describe also interparticle material exchange by 
simulating diffusive fluxes between particles by virtue of 
differences in chemical potentials, thus allowing the method 
to deal with a greater variety of systems. 

The algorithms of GenDPDE-M are expressed as Langevin 
equations for the relevant dynamic properties, which are, in 
this case, momentum, internal energy, and particle 
composition. 


Coupling between energy and mass exchange 


The original formulation of GenDPDE-M assumes that 
material exchange between mesoparticles is not coupled with 
energy exchange, which represents a limitation. In the present 
work, coupling between these two properties is introduced, 
taking into account that Detailed Balance not only fixes the 
amplitude of the random contributions, but also imposes 
Onsager’s Reciprocal Relations at the mesoscopic level. 
Equilibrium simulations employing the van der Waals 
Mixture model for a two-component system for the particle 
internal thermodynamics are used to validate the internal 
consistency of the stochastic model (Figure 1). Furthermore, 
simulations at non-equilibrium conditions are conducted to 
demonstrate that the model is capable of reproducing and 
controlling the Ludwig-Soret Effect, which depends precisely 
on the coupling between energy and matter transported at the 
mesoscale (Figure 2). 
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Figure 1: Equilibrium distributions of the particle temperature (top) 
and masses of the embedded species, Ar and Cl2 (bottom) (T=1.081, 
p=0.3152). The dotted and squared curves refer to coupled systems 
characterised by different values of the mesoscopic Onsager 
coefficients, whereas the continuous line refers to a decoupled 
configuration. Superimposition of the distributions indicates 
consistency of the proposed model. 
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Figure 2: Simulation box (top) and average profiles for the particle 
temperature (middle) and number concentration of one of the 
embedded species, Cl2 (bottom). The imposed temperature gradient 
along the x-axis between the hot (red) and cold (blue) slabs induces 
a concentration gradient in the system. The slope of these profiles is 
determined by the strength of the coupling between energy and mass 
transport, as shown by the differently coloured lines in the diagrams, 
which refer to configurations considering different values of the 
mesoscopic Onsager cross-coefficient. 


Conclusions 


In this work, we present the stochastic set of governing 
equations for the GenDPDE-M model when coupled energy 
and material exchanges occur, and show the consistency of 
the method together with several non-equilibrium results for 
the test system. With the addition of the coupling between 
energy and mass transfer, GenDPDE-M can be considered as 
the most general model for the simulation of complex fluids 
at the mesoscale when thermal fluctuations are relevant. 
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Introduction 


The Soret coefficient of electrolyte solutions features a non- 
monotonic dependence on the salt concentration. Early stud- 
ies of alkali halides reported a minimum below | M LiCl con- 
centration and 0 oC [J. Colombani, 1999]. The existence of 
minima has remained a matter of discussion for several years. 
More recently, non-equilibrium simulation studies reported 
minima in LiCl [di Lecce et al. 2017a, di Lecce et al 2017v], 
verifying previous experiments. Moreover, these simulations 
established correlations between the heat transport and the ob- 
servation of minima and the possible impact of these minima 
on thermoelectric effects. 


Very recently, a new set of experimental studies of iodide 
salts, amongst other electrolytes, provided convincing evi- 
dence for the minimum at a specific concentration, 1 M. A 
phenomenological explanation was proposed to rationalize 
these results [Mohanakumar et al., 2022]. 


The studies discussed above have focused so far on aqueous 
systems. However, non-monotonic dependencies might be 
expected in non-aqueous systems, and some experimental ev- 
idence supports this idea. Interestingly, minima in the Soret 
coefficient have recently been reported using non-equilibrium 
molecular dynamics simulations in combination with simple 
binary liquid mixtures modelled with the Lennard-Jones po- 
tential. That work shows convincing evidence for a minimum 
at near 50:50 composition. Furthermore, a correlation was 
proposed between the observation of the minimum and the 
non-ideality of the liquid mixture. 


This presentation will provide an overview of current studies 
of the non-monotonic dependence of the Soret coefficient 
with the composition of liquid mixtures and solutions and dis- 
cuss possible microscopic mechanisms responsible for this 
observation. 
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Introduction 

The DCMIX project is a milestone in the research of 
thermodiffusion phenomenon in multicomponent mixtures, 
and the key to develop an on-ground nucleus for diffusion and 
thermodiffusion research. Nevertheless, the lack of 
knowledge of the phenomenon in different types of mixtures 
raises the need for further analysis of the transport properties 
of binary mixtures coupled to ternary systems. The former are 
key to understand the fundamental mechanisms governing 
mass transport phenomena in multicomponent liquid 
mixtures, as the initial link to study several properties such as 
the density, viscosity, the shape and size of molecules etc., or 
the development and validation of both theoretical and 
numerical models. 

In this context, this work presents isothermal and non- 
isothermal transport properties as molecular diffusion, 
thermodiffusion and Soret in a hydrocarbons system of 
unequal molecular shape are analysed; the polycyclic 
aromatic hydrocarbon methylnaphthalene (MN), _ the 
monocyclic aromatic hydrocarbon toluene (Tol) and a chain 
shape alkane n-decane (nCio). Experiments were conducted 
in 15 binary and 10 ternary mixtures. Thus, an integral and 
transversal study of the system was carried out, characterising 
first and second order thermophysical properties from binary 
to ternary mixtures, covering the whole pool of concentration. 
In short, thermodiffusion, molecular diffusion and Soret 
coefficients of twenty-five binary and ternary mixtures of 
MN|Tol|zCio system at 25 °C are analysed together with all 
necessary thermophysical properties. 


Analysed systems 

Constituents of the system are arranged by the hydrodynamic 

principle on the descending order of density. Hereinafter, MN 

(c,), Tol (cz) and nCio (c3). For each pair of components, 

0.20, 0.40, 0.50, 0.60, 0.80 in mass fractions were analysed. 
0.0 


—- —- = 
0.0 0.2 0.4 0.6 0.8 1.0 
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Figure 1: Ternary diagram of MN|Tol|7Cio system. Light and 
dark green state points indicate the 10 ternary and 15 binary 
mixtures under study, while the blurred light blue refer to 
additional measurement points for the determination of 
thermophysical properties. Arrows determine the components 
direction within the triangle. 


For the ternary mixtures, in turn, we focused on the following 
concentrations: 0.80|0.10)0.10, 0.60/0.20|0.20, 0.40|0.40/0.20, 
0.40|0.20|0.40, 0.33|0.33|0.33, 0.20/0.60|0.20, 0.20|0.40|0.40, 
0.20|0.20|0.60, 0.10|0.10)0.80, 0.10|0.80)0.10. 


Experimental setups 

The thermophysical properties measured were the following: 
density, thermal and mass expansion coefficients, dynamic 
viscosity, thermal conductivity and thermal diffusivity. The 
former three were measured using the Anton Paar DMA 5000 
M density meter, changing either sample concentration or 
mean working temperature. On the other hand, the AMVn 
micro viscometer from Anton Paar was used to study the 
viscosity of mixtures, and the Thermtest THW-L1 instrument 
to measure thermal conductivity and thermal diffusivity, 
based on the transient hot-wire technique. The data on 
refractive index and derived solutal contrast factors were 
obtained by the Anton Paar Abbemat Multi-Wavelength 
refractometer. The instrument is operating at seven different 
wavelengths of visible light, 2 = 435.8-656.3 nm. 

The steady-state concentration separation in a 
thermogravitational column (TGC), together with 
aforementioned thermophysical properties, permits 
straightforward determination of the thermodiffusion 
coefficients. In this work, the extraction column and the 
thermogravitational micro-column were used. In the first one, 
samples are extracted from intakes at different column heights 
and the concentration is determined by density and/or 
refractive index measurements. In the second one, optical 
digital interferometry is used to track the constituents’ 
separation during the entire experiment. This method, allows 
determining the relaxation time of the mixture at the same 
time, and thus, the molecular diffusion coefficient in the case 
of binary mixtures. Molecular diffusion coefficients of the 
three symmetric binary-subsystems and ternary systems was 
determined by the Sliding Symmetric Tubes (SST) technique. 
Further information about the TGC and SST techniques can 
be found elsewhere (Larrafiaga et al. 2015, Seta et al. 2022). 
The Soret coefficients were determined by the indirect 
method, given by the following expressions for both binary 
and ternary mixtures, 


Sr = Dr /D, (1) 


(2) 


The analysis of a binary mixture in the TGC makes it possible 
to calculate the Soret coefficient by, combining the separation 
in the transient and steady state regimes. In the case of ternary 
mixtures, experiments in the SST were combined with the 
coefficients measured in the extraction column. 
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Results 

The thermodiffusion experiments in binary subsystems reveal 
positive thermodiffusion coefficient for all mixtures. In the 
case of the ternary mixture, MN moves towards the cold wall 
and is enriched in the bottom of the column. On the contrary, 
nC1o moves towards the hot wall and concentrates at the top. 
The thermodiffusion coefficient of the toluene, intermediate 
component, varies depending on the concentration analysed. 
Figure 2 shows the variation of concentration 0.20|0.20|0.60 
in mass fraction ternary mixture. 
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Figure 2: Variation of constituents’ concentration along the 
height of the TGC for MN|Tol|7Cio mixture at 0.20|0.20/0.60 
mass fraction at 25 °C. 


The transient separation and steady-state analysis (Figure 3) 
reveals the relaxation time and concentration variation of the 
mixtures, and therefore, the thermodiffusion, molecular 
diffusion and Soret coefficients of binary mixtures. 


0.04 


con = -1.459 z +0.518 
~ 0.48 Hy Soa | 
0.000 0.005 0.010 0.015 0.020 0.025 
Column Height (m) 
0 200 400 600 800 
Time (min) 

Figure 3: Result of a thermodiffusion experiment of MN|Tol 
0.50 gg" in the thermogravitational microcolumn. Relaxation 
time and separation in the steady state at 25 °C. The inset 
shows the change in the MN concentration with TGC height. 
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Figure 4: Concentration variation in the SST experiment for 
MN|Tolj7Cio mixture at 0.20|0.20|0.60 mass fraction at 25 °C. 


Tarragona (Spain), May 2023 


In the case of the ternary mixtures, the concentration variation 
of the two independent components under two different initial 
conditions allows the determination of the four coefficients 
that make up the diffusion matrix. Figure 4 shows an example 
of an SST experiment. 

Our experiments revealed that the Soret coefficients of MN 
are positive for all the analysed ternary mixtures while the 
ones corresponding to nCjo, are all negative. The Soret 
coefficients of toluene can be either positive or negative 
depending on the mixture. 


Conclusions 
The objective of this research is to determine the transport 
properties along the three binary borders, which allow 
generating knowledge about the ternary mixtures. In order to 
validate these predictions, the properties have to be measured 
experimentally. In the case of the thermodiffusion coefficient, 
Blanco et al. 2010 developed a combination rule, which 
predicted the thermodiffusion coefficients of ternary mixtures 
from binaries: 

Dri Lc Dy cic; LL DP cc. = 

am aij Qik 

where D7; is the thermodiffusion coefficient of component 
i of a ternary mixture, v,,. is the corresponding dynamic 
viscosity and a is the thermal expansion coefficient. D7’ 
and D4 are the thermodiffusion coefficients, Vij and Vix 
are the kinematic viscosities and aj; and aj the thermal 
expansion of the corresponding binary mixtures. c;, cj and 
Cx correspond to the mass fraction of each constituent in the 
ternary mixture. Thus, considering the number of mixtures 
measured around the system, the combination rule was tested 
and validated with a mean deviation of 15 % in this system. 
In the case of the Soret coefficients, following the work 
presented by Mialdun et al. 2021 the Soret vector model was 
used to verify the Soret coefficients of the ternary mixture. It 
was concluded that MN|Tol|nCio ternary system can be 
considered a nearly ideal mixture as THN|IBB|nC12. 
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Introduction 


Thermodiffusion, also called the Soret effect, is observed to 
influence all kinds of matter, from gases, liquids, to solids. In 
metals, the effect has been employed for nuclear enrichment 
and crystal growth, it is relevant for solder and manufacturing 
of integrated circuits, and has recently been shown to be of 
use in creating metallic nanowires. 

There exist several theoretical models attempting to predict 
binary thermodiffusion, but there is still no comprehensive 
model able to predict the Soret coefficient for a wider range 
of systems, as recently shown by (Hoang et al. 2022). The 
authors of such models often voice the need for more 
experimental data, which is very scarce in the case of metallic 
melts. The same holds for modeling thermodiffusion using 
computational simulations, where the dynamics of 
thermodiffusion in the simulations (dependence of 
thermodiffusion on mixture size ratio, dependence of Soret 
coefficient on concentration, etc.) should be compared with 
measurements. 

Almost all previous publications on liquid alloys are on 
binary systems where the atomic mass ratio of the two 
components is around two or less. The only exceptions so far 
are measurements on carbon and trace amounts in a solvent. 
The experimental data so far indicate that the Soret coefficient 
is generally dependent on the relative molecular weights of 
the species in the mixture. Measuring the thermodiffusion in 
liquid Al—Ag for different concentrations therefore provides 
a new insight into the dynamics of thermodiffusion in an 
atomic fluid with a high mass ratio (Srinivasan et al. 2013). 
Here, we study thermodiffusion in the liquid alloy Al-Ag as 
a function of composition using a newly developed 
experimental setup. 


Method 


Using X-ray radiography (XRR), we have measured the 
thermodiffusion in liquid Al-Ag, with compositions ranging 
from AlgoAg2o to AlsoAgso. The samples are around 10 mm 
long and have a diameter of 1.3 mm in the liquid state. 

The setup that is used supports four samples, of which two 
samples are reference samples to calibrate the relation 
between image brightness and sample concentration. The two 
reference samples are chosen at 5 percentage points higher 
and lower than the main samples. This also means that 
reference sample data for each X-ray image are available, 
which are beneficial in case of fluctuations in the X-ray 
source over the time of the experiment. 

The samples were melted by heating to 1023 K. After 
homogenization for about two hours in an isothermal state, a 
temperature gradient of around 10 K/cm was established by 
switching off the heater at the bottom of the samples. 


Results and discussion 


Figure | shows the concentration difference for an AljsAgos 
sample over time with the hot end fixed at 1023 K. The figure, 
displaying the concentration difference of silver between the 
hot end and the cold end of the sample, shows that the silver 
concentration is reduced in the hot end, and increasing at the 
cold end, meaning that silver migrates to the cold side. 
Using the concentration gradient data from the non- 
isothermal equilibrium phase (highlighted on the right in the 
aforementioned figure) together with the last data of the 
isothermal phase (highlighted on the left in the 
aforementioned figure), the Soret coefficient can be 
calculated directly using 
ACq, 
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. AT co(1 — Co) 
The Soret coefficient is determined to be (0.940.3)x10% K"! 
(Kriger et al. 2023). 
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Figure 1: Concentration difference across a sample of Al7sAg2s. The 
vertical line indicates when the temperature gradient was initiated. 
The exponential curve is the fit of the equation for the time- 
dependent concentration separation to the experimental data after the 
temperature gradient had stabilized. The colors indicate the basis for 
the averaged data used in calculation of the Soret coefficient. 


Thanks to the time-resolved information from in situ XRR, 
the interdiffusion coefficient can simultaneously be measured 
(Sondermann et al. 2019). By looking at the total 
concentration difference Ac over time, its time-dependence 
can be fitted to the data in figure 1, where the equation is 
given by 

8 het e7(2nt1)7t/6 

bel =Aee (1-5 ae 
rl 

with the characteristic time 0 = — for a sample of length 
L and of a material with diffusion coefficient D. We then get 
a diffusion coefficient of (3.8+1.5)x10° m?/s. This is in good 
agreement with previous measurements by (Engelhardt et al. 
2016) on the diffusion coefficient for the Al-Ag system. 
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The similarity in the diffusion coefficients corroborates the 
measured Soret coefficient, as the observed process goes at 
the speed predicted by theory. 


Using thermophysical values from the literature, we can 
calculate the predicted Soret coefficient by the current 
theoretical models (Eslamian et al. 2010, Shukla et al. 1998). 
The predicted Soret coefficients of the two models are 
presented in figure 2. 


20 30 40 50 


Ag (at.%) 


Figure 2: Soret coefficients accumulated from all experiments, 
measured at a mean temperature of 1018 K and with temperature 
differences of around 10 K. Error bars show one standard deviation 
for the Soret coefficient for each concentration. The dashed line 
shows the weighted mean of all our measurements. The predicted 
values by the models of (Eslamian et al. 2010) and (Shukla et al. 
1998) are shown in orange triangles and green pentagons 
respectively. 


Both models predicted the same sign for the Soret coefficient 
of liquid Al-Ag alloys as found in our experiments, i.e., 
correctly predicted the direction of thermodiffusion for the 
two components in the alloy. But the numerical values were 
off by at least a factor of two with respect to our measured 
values, and also the strong concentration dependence 
predicted by the models was not observed. This shows that 
there is still need for a theory to describe thermodiffusion 
even in binary atomic liquids. 
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Conclusions 


Using X-ray radiography (XRR), the Soret and interdiffusion 
coefficients can be measured in binary liquid alloys. It also 
allows several samples to be analyzed simultaneously as well 
as several times in succession, with relative ease. We 
investiaged the liquid Al-Ag system with compositions 
ranging from AlgopAg20 to AlsoAgso. A Soret coefficient of 
(0.9+0.3)x10° K-! was found over the entire measured 
composition range, with silver migration to the cold side, 
which is of the same order of magnitude as previous 
measurements on similar binary systems. Using experimental 
data from the literature, models to predict the Soret coefficient 
were tested. The closest model is off by a factor of two with 
respect to our measured values. The method of XRR for in 
situ measurements at different times also allows the 
measurement of the interdiffusion coefficient, which was 
found to be in accordance with previous interdiffusion 
measurements on the Al-Ag system. In the future, we want to 
analyse thermodiffusion in liquid Al-—.In, which has a 
miscibility gap, where the derivative of the chemical potential 
with respect to concentration (a crucial parameter in models 
for thermodiffusion) is sensitive to concentration/temperature 
changes. Ternary alloys will also be investigated. 
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Introduction 


The geothermal gradient within deep saline aquifers used for 
carbon sequestration may result in an uneven spatial 
distribution of CO2 molecules. Nonetheless, the investigation 
of thermodiffusion of CO2 in aqueous mixtures remains to be 
thoroughly pursued The first attempt to quantify the 
thermodiffusion in CO2-H2O mixture has been made by 
Windisch et al. (2012). The authors applied Raman 
spectroscopy to quantify the concentration gradient of CO2 
due to a temperature gradient within a temperature range 
between 288 and 318 K and at pressures to 10 MPa. They did 
not observe detectable Soret effect, within the uncertainty of 
the measurements. Guo et al. (2018) have also performed 
experiments using UV Raman spectroscopy in a temperature 
gradient established between 295 and 353 K and pressures of 
20 and 30 MPa. The authors have examined a saturated 
mixture of CO2z, and the observed concentration gradient 
established was probably due to the temperature dependency 
on solubility rather than the thermodiffusion. 


Thermodiffusion can also be investigated by molecular 
dynamics (MD) techniques (Galliero et al. 2008, Nieto- 
Draghi et al. 2005, Perronace et al. 2002). With Boundary- 
Drive Non-Equilibrium Molecular Dynamic (BDNEMD) 
simulations, heat flux or a temperature gradient is generated 
in the simulation box, and, by the system response, a 
concentration profile is established. The thermal diffusion 
factor a7; can be computed to quantify the phenomenon. For 
a binary system, az; is defined as: 

r T dx; 1 
x;(1 = x;) aT ( ) 
where T is the absolute temperature, x; is the composition in 
mole fraction of component i, and J is the thermodynamic 
factor that should be accounted for in non-ideal liquids I” # 
1) because diffusion is driven by chemical potential gradient. 
For multicomponent mixtures, the thermal diffusion factor 
definition is not straightforward, and there is dependency on 

the transport diffusion coefficients. 


ari = 


In this work, we evaluate the thermodiffusion of CO: in saline 
solutions using non-equilibrium molecular dynamics 
simulations. First, the Soret effect for CO2-H2O mixture is 
computed. Then, salts are added to the system, and the 
thermodiffusion in the ternary mixture is investigated. 


Methods 


The eHeX method is the BNEMD method chosen 
(Wirnsberger et al. 2015). The heat flux is generated by 
scaling up and down the velocities of particles that are 
located, respectively, in the hot and cold predefined regions 
of the simulation box as illustrated by Fig. 1. 
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Figure 1: Representation of the simulation box with cold 
and hot regions. 


Different combinations of force fields are selected to 
investigate the sensitivity of the results. Water molecules are 
represented by the SPC/E and TIP4P/2005 models, whereas 
COz by EPM2 and TraPPE. The simulation box is set with 
5000 molecules, and a density of p = 1000 kg.m*?. The CO 
composition in mole fraction is set at 0.02; a value below the 
solubility at the specified pressure and temperature. The mean 
temperature and the temperature difference between the cold 
and hot regions are 350 K, and 50 K, respectively. All 
simulations are performed in LAMMPS. First, an NVT 
equilibration is performed for 1 ns to keep the system close to 
the desired mean temperature, and then, the eHeX method is 
turned on for 80 ns. The thermodynamic factor is computed 
from equilibrium molecular dynamics simulations through 
the Kirkwood-Buff integrals using the method from Dawass 
et al. (2019). The saline effect on COz2 thermodiffusion is 
investigated by adding NaCl in the mixture at a composition 
of 2 molal, and then evaluating the segregation of the 
components in the presence of the temperature gradient. 


Results 


Under the analyzed conditions, for the binary CO2-H2O 
mixture, CO2 accumulates in the cold region (arco2 > 0), and 
the thermal diffusion factor is around 1 (Coelho et al. 2023). 
The same result is reproduced by all force field combinations; 
hence only the conventional SPC/E + EPM2 combination will 


be presented (Figure 2). 
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Figure 2: Temperature and CO2 composition profile in the 
simulation box at T= 350 K, p = 1000 kg.m®, xco2 = 0.02. 
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Water molecules interact through highly directional hydrogen 
bonds. The hydrogen bond network is disturbed by increasing 
the temperature. CO. molecules do not have a net dipole 
moment, but they interact through van der Waals and 
quadrupole-quadrupole interactions. For CO2, the higher the 
temperature, the higher are these interactions. The cross- 
interactions between CO2-H2O molecules could be through 
tetrel bonds or hydrogen bonds. From the radial distribution 
functions, the tetrel bonds are more significant than the weak 
hydrogen bonds for this mixture and are decreased by 
temperature, whereas the hydrogen bonds are enhanced. 


The arco2 dependency on temperature is evaluated by 
performing simulations at various mean temperatures at the 
same pressure as the initial condition (400 bar). The thermal 
diffusion factor decreases by almost 50% by increasing the 
temperature from 300 to 350 K, and is almost negligible at 
400 K, where there may be a sign change (Figure 3). The 
pressure effect is evaluated by equilibrating the mixture at 200 
bar; however, the pressure influence on the Soret effect is not 
as pronounced as the temperature (arco2 decreases by 28% by 
decreasing the pressure from 400 to 200 bar). 


@r,co, = -0.017T + 6.9 


300 320 340 360 380 400 
Temperature / K 


Figure 3: Thermal diffusion factor of CO2-H20 vs. 
temperature at P = 400 bar, and xco2 = 0.02. The dashed line 
represents the linear interpolation. 


Overall, the separation between CO2 and water molecules due 
to a temperature gradient is enhanced by increasing the 
mixture density (lower temperature and higher pressure). 
CO, accumulates in the cold region where its interactions 
with other CO, molecules are reduced, whereas water 
accumulates in the hot region where it has reduced self- 
association. This indicates a CO 2 accumulation in the top 
region of the aquifers. 


In brine mixtures with CO2, NaCl salt accumulates in the cold 
region. The repulsive interactions between CO2-NaCl leads to 
a less thermophobic behavior of CO2. At 350 K, there is 
almost no detectable Soret effect, and at 400 K, CQ2 
accumulates in the hot region. 


Conclusions 
The focus of our work is the investigation of the 


thermodiffusion of the CO2-H20 mixture at subsurface 
conditions. Under a temperature gradient, CO2 accumulates 
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in the cold region, where the hydrogen bond network is more 
established. The Soret effect for CO2 in water is enhanced as 
density increases by decreasing the temperature or increasing 
the pressure. The temperature effect on the thermal diffusion 
factor is pronounced. The electrolytes have a pronounced 
influence on CO: distribution, changing the thermal diffusion 
factor sign at a lower temperature. 
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Introduction 


Introduction 

Most mixtures in nature and industry are commonly 
multicomponent and subject to non-equilibrtum conditions, 
thus the knowledge of the transport phenomena occurring in 
them is of great interest for many applications. When a liquid 
mixture is subject to a temperature gradient, its constituent 
species move, and a concentration gradient is generated. The 
complex interaction between heat and mass transfer processes 
can result in a variety of dynamical behaviors due to 
instability, being the coupling between temperature and 
concentration gradients, also known as thermodiffusion or the 
Ludwig-Soret effect (Soret, 1879), one of the mechanins that 
can drive convective instability. Soret effect is quantified by 
the Soret coefficient Sy, the ratio between the thermodiffusion 
coefficient Dr and the molecular diffusion coefficient D 
(Kohler et al., 2016). Past years brought a great development 
on optical methods for the measurements of both molecular 
and thermal diffusion coefficients for binary mixtures 
(Mialdun et al, 2011) motivated by their high accuracy allied 
to anon distressing feature on the diffusive process. However, 
the fact that the sign of the Soret coefficients of the various 
species in a mixture can be different, disturbing the system, 
makes the study thermodiffusion in complex mixtures more 
challenging. Consequently, only a limited number of 
thermodiffusion coefficients for binary and ternary mixtures 
are available in literature and, often, data are not consensual. 
In this work, we examined the Soret coefficient of the mixture 
composed of a nonpolar hydrocarbon solute triethylene glycol 
(TEG) and water, exploring the full concentration region of 
the mixture at 25 °C. Experimental method for the 
determination of Soret coefficients is based on the optical 
digital interferometry (ODI) technique (Mialdun et al, 2011). 
Research on triethylene glycol is particularly important due 
to its applications as a dehydration agent in natural gas 
streams to prevent corrosion in pipelines (Trueba et al., 2022). 
Moreover, this is a non-ideal mixture, and thus, it has a 
potential to demonstrate a complex dependence of the Soret 
coefficient on concentration, including its sign change. 


Experimental 


A. Materials 

The reagent grade triethylene glycol (purity 0.99%, CAS 
Number: 112-27-6) and water extra pure, deionized (CAS 
Number: 7732-18-5) were obtained from Acros Organics and 
were used as received, with no further purification. Liquid 
samples of ~15 g were prepared in mass fraction, by weighing 
each component using the Sartorius 1712 analytical balance 
with a resolution/capacity of 0.1 mg/160 g and then remixed 
by a magnetic stirrer during several hours. 


B. Experimental 

interferometry 

The experimental setup for the measurement of the Soret 

effect used in this work was optical digital interferometry 

(ODD, a validated and reliable optical method very well 

described in literature (Mialdun et al, 2011; Santos et al. 
2022). 


setup: Optical digital 


Results 


A. Overview of the component separation 

The study of the separation in the TEG—water mixtures was 
done over the complete composition range of TEG, using a 
0.1 mass fraction step. In the low TEG concentration range, 
the mixtures featured stable separation, well accessible by 
ODI, with the maximum optical signal found at c ~ 0.2 mass 
fraction. A further increase of the TEG concentration resulted 
in a gradual decrease of the optical signal corresponding to 
steady-state separation, which was found unstable (negative) 
in the experiment at c = 0.5 mass fraction. Several 
experiments, made at even higher concentrations of TEG (at 
c = 0.75 and 0.98 mass fraction), confirmed the instability of 
the separation over all regions of high TEG concentrations, 
until c = 1. The summary plot of the observed concentration 
dependency of the optical separation is shown in Fig. 1. 
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Figure 1: General trend of separation for the mixture TEG—water 
(dots represent the experimental data, and the dashed line is the fit). 


Fitting of the results with a first-order R-K polinomia, which 
considers the disappearance of the separation at the diluted 
limits, (Redlich et al, 1948) provides a reasonable description 
of the explored region and a satisfactory guess about the 
unexplored one. The obtained fitting curve, represented by the 
dash line in Fig. 1, allowed to precisely locate the 
concentration specific to the sign change of the separation, at 
a 0.45 mass fraction of TEG. Until this point, the separation 
of components is stable and develops normally, having the 
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point of maximum separation at 0.2 mass fraction of TEG. In 
an attempt to improve the accuracy of the dependence of 
separation, new measurements at intermediate concentrations 
in the region of stable separation, with a smaller step of 0.05 
mass fraction, were added. 

While the ODI setup configuration is optimized for the 
measurement of mixtures with the positive Soret separation it 
was still possible to measure two extra data points that belong 
to the region with unstable separation, and the extraction of 
data is explained with detail in next point. 


B. Unstable separation 

When running the experiments with mixtures featuring small 
unstable separation (e.g., at c = 0.50 and 0.98 mass fraction), 
we made an important observation. At first, the separation, 
even being generally unstable, develops in accordance with 
the analytical solution. At later time, when enough amount of 
the heavier component (TEG) has accumulated at the top of 
the cell, the 1-D character of the developing separation breaks, 
and a hydrodynamic instability develops (Santos et al., 2022). 
The instability should start after a certain threshold and may 
take either long-wave or finger-type character, with both the 
threshold and the character being dependent on 
thermophysical and transport properties of the mixture, as 
well as on the applied temperature difference and the cell size. 
The most important feature here was the duration of the step 
of the stable development of the separation. In the two most 
interesting cases, at c = 0.5 and 0.98 mass fraction, this 
duration was found as 2—3 h. The separation development in 
the case of c = 0.5 is shown in Fig. 2. 


separation 
develops 
normally 


an/az / 10-3 m-} 
i) 
°o 


. 
0.0 4 Py meee Toten, . 
7 1 . . eee 
! 
4 t 
0 


i 


an a A RT 
Figure 2: General trend of separation for the mixture TEG—water of 
c = 0.5 (dots are experimental data, and the dashed line is the fit). 
The wrapped optical phase map shown below corresponds to the 
beginning of instability development. 


It is seen in the plot that, during the first hours, the separation 
develops in a perfect agreement with the 1-D analytical 
solution and that due to the faster image acquisition at the 
beginning of the experiment, this time interval is quite 
densely covered by interferograms. 

It appeared that this amount of data is enough for fitting and 
extracting the separation magnitude and diffusion coefficient, 
despite that the separation does not reach even a half of the 
steady-state level. The 2-D phase map specific to the end of 
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the short period of the stable separation development is shown 
in Fig. 2, below the separation curve. There are visible traces 
of the developing instability in the corners of the cell, while 
in the liquid bulk in the center of the cell, the distribution 
resembles the steady one (Santos et al., 2022). Thus, in certain 
cases of unstable separation, the data allow access to the 
transport properties in the same manner as for fully stable 
separation. This, however, is not always the case. For 
example, the experiment with the mixture of c = 0.75 mass 
fraction concentration did not allow one to apply this 
approach since the instability development started after 0.5 h, 
obstructing a further extraction of meaningful data. 
Nevertheless, this possibility let us add two additional points 
to our results, having a reliability comparable with all other 
cases featuring the stable separation. 


Conclusions 


Soret coefficients St were determined for the TEG—water 
mixture with ODI setup, showing that St changes its sign at 
0.1 mol fraction of TEG (0.45 mass fraction of TEG) and that 
the separation of the components at high concentration of 
TEG is unstable. 
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Introduction 


In the gravity field or under microgravity, pure thermo- 
diffusion leads to very weak species separation. To increase 
the species separation in the presence of gravity, many authors 
use thermo-gravitational diffusion in vertical or inclined 
columns (TGC), (Platten et al. 2003), (Charrier-Mojtabi et al. 
2007) and (Seta et al. 2019). The flow velocity strongly 
depends on the temperature difference, AT, imposed between 
the two walls facing each other, and also on the thickness, H, 
between them. Indeed, an optimum species separation cannot 
be obtained in relation to AT and H simultaneously. This 
problem admits an optimum which is only function of H for 
fixed AT or function of AT for fixed H, (Mojtabi et al. 2019). 
The coupling between shear-driven convection and thermo- 
diffusion is a complex phenomenon due to the interactions 
between the different forces applied to the fluid mixture. 

Up to now, the species separation obtained in vertical columns 
has been limited. The convective flow velocity depends on the 
temperature difference AT, maintained on the vertical 
surfaces of the column and the thickness H separating these 
two surfaces. For a given thickness, to increase the 
importance of thermodiffusion the temperature difference AT 
imposed on the two vertical surfaces must be increased. This 
results in an increase in the convective velocity of the fluid 
resulting in a decrease in species separation. 

For a fixed AT value, the optimum species separation is 
obtained for a vertical column thickness less than one 
millimeter. 

In this work, the authors sought to dissociate the fluid flow 
allowing the species separation from the temperature 
difference imposed on the two parallel surfaces of the cavity. 
In addition, appropriate thermal boundary conditions were 
determinated in order to increase the temperature difference 
without generating natural convection for a positive or 
negative Soret coefficient of the binary fluid. 

It was shown that it is possible to significantly increase the 
importance of the species separation for a binary mixture, 
using forced convection obtained with a uniform translation 
displacement of the horizontal isothermal walls maintained at 
Trop and Tyot respectively, with AT = Trop -Tyot- 

We verified that it is possible to obtain a significant species 
separation whatever the sign of the thermodiffusion 
coefficient of the binary fluid and this procedure can be 
extended to multi-constituant fluids. In addition, this 
procedure can be used with cell thickness greater than the 
vertical column thickness leading to optimal separation. 

This procedure would also make it possible to obtain 
measurements of negative thermodiffusion coefficients, in the 
earth laboratory, without using microgravity experimentation 
( Schram et al. 2021) 


Figurel: Schematic of the rectangular cavity 


1. Analytical solution in the case of shallow cavity 

We considered a rectangular cavity of large aspect ratio 
A=L/H>>1, where H is the height of the cavity along the y 
axis and L is the length along the x axis. The cavity was filled 
with a binary mixture of density p and dynamic viscosity i. 
The corresponding dimensional boundary conditions are: 


pit aw € [0, H] 
v= UL: aa "ax ax , y , 
y=0,H 
V > = Oc _ OT 
V = fUpéx, Upéx,T = Toot» Teop:D 3 + Dr zy = 0 


where D7 = Cy(1 —Co)Dr , Cy and Dr are respectively the 
reference mass fraction and the thermodiffusion coefficient. 
In the case of a shallow cavity A >> 1, the parallel flow 
approximation, used by many previous authors, was 
considered. The streamlines are all parallel to the x axis 
throughout the cavity except for the vicinity of the insulated 
walls x =0 and x = L, which gives: 


V, = U,()éx Tp = bx + g(y) and C, 
= mx + f(y) 


The constants b andm _ represent respectively the 
temperature and mass fraction gradients along the x axis. 
The constant b is zero since constant temperatures are 
imposed on the walls y = 0,H. Replacing in the mathematic 
formultation of the problem, the velocity, the temperature and 
the mass fractions by their expression U,,T, and C, , we 
obtained the analytical solution of the simplified equations. 
The mass fraction gradient was obtained by writing that the 
mass flux associated with the constituent of mass fraction 
C, through any cross-section of the rectangular cavity 
perpendicular to the x-axis is equal to zero: 
" (Uply — Do = C(1 — &) 2) ay = 0 
[ate ~ D5? — cu - 6) Say = 

The extremum of m, with respect to U;, and the coefficient f 
is obtained for f = —1. U, and mare the solutions of fourth 
degree algebraic equations leading to two real solutions and 
two complex solutions for the binary mixture studied. The 
two real solutions are opposed (Uj, -m) and (-U,, m) and lead 
to the same physical solution. 

We also conducted a numerical study of this problem, without 
using any restrictive hypothesis. A comparison between 
analytical and numerical solutions was performed for two 
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binary mixtures, one with D; > Oand the other with Dr < 0. 
These comparisons showed a good agreement between the 
analytical and numerical results. 


2. Results obtained for two binary mixtures with 

D; >Oand D; <0 

How to obtain forced convection alone, within a horizontal 
binary fluid layer, in the presence of the gravity field? A 
distinction should be made between binary fluids with a 
positive thermodiffusion coefficient Dr and negative Dr. 


a- Ethanol water mixture with 60.88% water, with Dr > 0 
For the Water-Ethanol mixture, the densest constituent, 
namely water, has a Dr; >0. In this case, it would be 
necessary to operate with a layer of fluid heated from above. 
Indeed the densest constituent migrates towards the cold 
bottom wall, the binary fluid layer is thermally infinitely 
stable. As U, values leading to optimal separations are of the 
order of 10 > m/s, forced convection flow will remain stable 
even for U, values well above 10 ~ m/s. 

The displacement along the x-axis of the walls maintained at 
Tpot and T;o,, combined with the thermodiffusion leads to 
species separation between the two ends of the cell (x = 
0 and x = L), for AT = Troy — Tyoe > 0 and for H2Jmm. 


0.63 
0.625 
0.62 
0.615 
0.61 
0.605 
0.6 
0.595 
0.59 


Figure 1: Mass fraction field, C,, for Water-Ethanol binary mixture. 


We present some results obtained for Water-Ethanol, 60.88% 
water for H = 2.10°?m,L=0.5m and AT = 10°C: 
Numerical: U, = +3.2107-5m/s; m= —0.521 m+ 
Analytical: U, = +3.2 10°°m/s; m = —0.514m7? 

When increasing AT, the velocity U, and the mass fraction 
gradient C, increase. When H increases at fixed value of 
AT, the mass fraction gradient decreases. 


b- Water-Isopropanol mixture with 90% water with Dr < 0 
For the Water-Isopropanol mixture, the densest constituent, 
namely water, has a D; < 0. In this case, this constituent 
migrates to the warmer wall. Then it would be necessary to 
operate with a layer heated from below, in order to prevent 
the onset of natural convection. This corresponds to the 
thermal configuration of Rayleigh-Bénard, but here the 
densest constituent is more concentrated at the bottom of the 
horizontal layer. Natural thermal convection can only occur 
for thermal Rayleigh values greater than the linear thermal 
critical Rayleigh number, Ra, = 1708. Natural convection, 
for the Rayleigh-Bénard configuration in the presence of 
binary mixture, can only occur for temperature differences 

a : 0 i 8 

Q ie) ie) is) i) 


Figure 3: Mass fraction field, C,, for Water-Isopropanol binary 
mixture. 


AT =Thot — Trop > AT, associated with Ra, = 1708, 
(Charrier-Mojtabi et al. 2007). 
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We present some results obtained for water-isopropanol 
mixture with 90% water for H=1.107-?m and AT = 
10°C: 

Numerical: U, = 7.81.107° = ; m=—0.462m71 
Analytical: U, = 7.81.10-°=; m = —0.467 m™ 

When increasing AT, the velocity U;, and the mass fraction 
gradient, C;, increase. 

For H =2.10-3m and according to the values of the 
thermophysical characteristics of binary mixture studied, the 
thermal Rayleigh number could become much higher than 
Ra, = 1708 leading to the onset of mixed convection. 


Conclusions 


This procedure of species separation using forced convection 
in the presence of the gravity field allows us to: 
-considerably improve the species separation (several orders 
of magnitude greater than the one obtained in vertical 
columns, TGC). 

- to be able to measure, in the Earth laboratory, not only 
diffusion and thermodiffusion coefficients for mixtures with 
positive thermodiffusion coefficient but also for mixtures 
with negative Soret coefficients. The measurement of the 
Soret coefficient for these mixtures is extremely difficult to 
obtain in the earth laboratory. This led several teams within 
the framework of the ESA DCMIX project to measure these 
coefficients in the ISS. 
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Introduction 


Fluid flow in macroscopic systems is always driven by 
pressure gradients or by external forces and is accurately 
described by continuum approaches, like the Navier-Stokes 
equation. However, starting from the investigations led by 
Derjaguin, it became clear that in sub-micron systems 
interfacial effects become relevant and additional 
mechanisms leading to a mass flow are possible (Derjaguin 
1941). 

A remarkable example is provided by thermo-osmosis, a 
surface-induced phenomenon, where the external “field” 
driving the fluid flow is a temperature gradient. 

Thermo-osmotic flows are expected to be minute on a 
macroscopic scale because the bulk fluid is set into motion 
due to an interfacial effect. However, the mass transport 
arising from any surface phenomenon becomes relevant every 
time one of the system’s dimensions becomes comparable to 
the correlation length or to the mean free path of the fluid. 
This situation occurs in several systems of technological and 
biological interest (Barragan 2017): Industrial miniaturization 
has made available a large variety of nano-channels, paving 
the way to the development ofnano-fluidics and many natural 
systems, such as membranes and gels, are characterized by 
porous networks with sub-micron diameters, where thermo- 
osmosis play an important role. 


State of the art 


Thermo-osmotic effects have been known since over a 
century and were investigated in rarefied gases starting from 
the studies by Maxwell (Maxwell 1879) and Reynolds 
(Reynolds 1879). The microscopic understanding ofthis out- 
of-equilibrium phenomenon in gases, first formulated by 
Maxwell within kinetic theories, is quite subtle: The flow, 
named thermal creep within this context, originates from the 
tangential stress on the gas layer near the confining surface 
due to the applied temperature gradient and is directed from 
the cold to the hot side. The effect is strongly related to the 
properties of the scattering event ofa gas molecule on the wall 
and disappears for an ideal gas confined by purely reflective 
hard walls. 

Conversely, in the liquid regime, the phenomenological 
descriptions usually adopted rely on macroscopic approaches, 
such as non-equilibrium thermodynamics and Navier-Stokes 
equations (Derjaguin 1941). However, the use of continuum 
theories can be justified when the relevant physical quantities 
vary on a length scale much larger than the typical range of 
the interaction: In the presence of a confining surface this 
condition is no longer satisfied, because, near the interface, 
the fluid properties eventually driving the phenomenon may 
display strong but short-ranged, modulations. Therefore, it is 
not surprising that the classical macroscopic paradigm can fail 


to predict even the direction of the induced flows (Piazza et 
al. 2008). 


Towards a microscopic theory of thermo-osmosis 


Recently, the interest on thermo-osmosis has been revived 
by molecular dynamics simulations aimed at a microscopic 
understanding of this effect (Wold et al. 1999, Galliero et al. 
2002, Fu et al 2017, Ganti et al. 2018). 

Stimulated by these findigs, the authors proposed a first- 
principle and unitary theory of thermo-osmosis (Anzini et al 
2019, Anzini et al. 2022). This approach, valid both for gases 
and for liquids, is based on a generalization of the linear 
response theory formalism to the case of an inhomogeneous 
fluid close to a surface. This theory sheds light on the 
microscopic mechanisms driving thermo-osmosis and 
quantitatively relates the fluid flow to the properties of the 
fluid-wall interface via suitable dynamical correlation 
functions at equilibrium. 

In this contribution we present a general analytical 
solution of these equations, which expresses the mass flow 
induced by a thermal gradient in terms of the mass current- 
heat current correlation function. This solution allows to 
confirm the interpretation of thermo-osmosis as an effect of 
the variation of the local enthalpy density induced by the 
presence of the wall, making contact with previous 
phenomenological approaches. 

To validate the theoretical predictions in a controlled 
environment, we performed extensive non equilibrium 
molecular dynamics simulations for a simple model of one 
component fluid and few wall-particle interactions. In most 
of the examples we will present, the walls are “passive”, i.e., 
they just act as an external force on the fluid molecules, 
orthogonal to the surface, thereby conserving the particle 
momentum parallel to the surface ina scattering process. 

The same microscopic theory allows to derive an 
analytical form of the velocity profile in the low-density 
regime when the confining surfaces violate the momentum 
conservation during the wall-molecule scattering. In this case, 
the comparison between numerical and analytical results is 
encouraging. 


Conclusions 


Thermo-osmosis is an interesting effect per se, being the 
simplest example of thermal force, and is expected to play a 
relevant role in different physical frameworks, from 
engineering to biophysics. 

In addition, one of the most important effects of thermo- 
osmosis probably occurs at the surface of colloidal particles 
immersed in a liquid or a gas, where the ensuing fluid flow 
pushes the colloidal particles through the fluid, giving rise to 
thermophoresis. 
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A microscopic study of thermo-osmosis is also 
instrumental for defining the correct boundary conditions (i.e. 
the slip velocity) for effective macroscopic approaches, based 


on hydrodynamics and the Navier-Stokes equation, 
describing fluid flow in confined systems. 
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Polymer solutions can be driven by temperature-dependent attractive interactions to organize into highly 
porous and tenuous gels displaying a very low elastic modulus. In this work, we investigating the optothermal 
effects generated by the weak adsorption of an IR laser beam on Mebiol, a block-copolymer hydrogel whose 
strength and elastic modulus grow with T that is extensively used as a substrate for cell growth. Structural 
changes induced the induced spatially-localized thermal gradient are detected by Photon Correlation Imaging 
(PCI), an optical technique that provides space-resolved information on the microscopic dynamics, allowing 
one to detect at the same time the presence of displacement and strain field. A schematic view of our PCI 
apparatus is shown in the figure. Preliminary measurements show that, even for moderate heating, a 
dramatic change of the microscopic gel dynamics, accompanied by the evidence of consistent strain effects, 
takes place. On switching off the laser, the gel slowly reverted to its original state, but only if the IR irradiation 
is kept for a limited time: For longer exposure, structural deformation became permanent, hence the gel 


displayed a plastic behavior. 
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Photon Correlation Imaging setup, including an IR laser beam 


operating at 980 nm for local sample heating 
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Introduction 


The thermal Marangoni effect, also called thermo-capillary 
convection, is fluid creep flow along a surface caused by a 
temperature gradient. The actual driving force is the gradient 
in surface tension. When a temperature gradient is applied to 
a fluid mixture in a porous medium, thermo-capillary 
convection and thermodiffusion may occur simultaneously. 
We investigate how the two processes couple and why 
thermo-capillary convection must be considered in 
measurements of Soret effects in porous media. 


If the permeability of the porous medium is low and the 
medium has sufficient mechanical strength, thermo-osmosis 
can generate a pressure gradient in the fluid. This pressure 
gradient can amount to several bars per degree [1]. We show 
how thermo-osmosis can be included in the theory for 
thermodiffusion. 


Combinations of type of fluid, type of porous medium, and 
effects of a temperature gradient are shown in Table 1. 


Table 1. Effects driven by a temperature gradient in different 
systems. 


Mixture 


nud One component 
Medium P 


Thermo-osmosis 


Thermo-osmosis Thermodiffusion 


Low permeability 


High permeability, 


bulk fluid Thermodiffusion 


No effect 


Theory 


The flux equations for a binary fluid mixture in a porous 
medium may be expressed as [2,3] 


; 1 1 xy 

Iq = Lag¥ (=) — Lav VP ~ Lap Vc (Ja) 
1 1 x4 

Iv = byg¥ (=) — Ly VP - Ly Vine (Ib) 
1 1 x4 

Jp => Log (=) Lov VP Lypp T Vic (1c) 


where Jj is the measurable heat flux and Jy and Jp are the 
volume flux and diffusion flux, respectively, defined as 


Jv =ShV1 + JoV2 (2a) 
lh te 
Jp = ee es (2b) 


were J; is the mass flux, V; the partial molar volume, and 
x; the mole fraction of component i , respectively. 


Furthermore, the L-coefficients are the Onsager coefficients, 
P is pressure, T is temperature, and yp, is the 
compositional contribution to the chemical potential from 
component 1. 


For a binary fluid mixture in a porous medium, we define the 
Soret coefficient as 
1 Vx. 

s= (Seq) 
X4xX2 VT 


(3a) 


Jv=Jp=0 


= 1 Lovlyq — LywvL pq (1 4 dlny, 
X4X2RT* \ LypLyy — Loy 


=4 
dOln = (3b) 


where R is the gas constant and y, is the activity coefficient 
of component 1. The condition Jy = Jp = 0 means no net 
flow through any cross section of the porous medium. This is 
achieved in molecular simulations by restricting the flow 
through the ends of the simulation box to keep a stationary 
state. 


The thermo-osmotic coefficient quantifies the thermo- 
osmotic effect. It is defined as 


a 
VI) jy=Jp=0 — \AT/ y=Jp=0 
1 (LovLpg — LopLvaq 
7 r( LypLyy — Lay ee 


where the differences AP and AT between properties in the 
bulk phases are used as an approximation. 


NEMD Simulations 


We will report results for thermodiffusion and thermo- 
osmosis in slit pores. Examples of the system layout for a one- 
component fluid of Lennard-Jones/spline particles is shown 
in Figure 1 (top). Non-equilibrium molecular dynamics 
(NEMD) simulations were done with the LAMMPS software 
[4]. Temperature gradients were established by heating the 
parts of the “hot fluid” and cooling parts of the “cold fluid”. 
The particle momenta of the left and right boundaries of the 
MD box were set to zero in order to ensure zero net mass 
fluxes in x-direction. The temperature gradients were the only 
external forces on the system. 


Preliminary results 
Result for the case shown in Figure | include the velocity 
profiles in the bottom panel of the figure. The wide pore has 


a distinct flow pattern with creep flow along the walls in the 
direction of high temperature. A consequence of the zero net 
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mass flow is that the fluid is forced in the opposite direction 
in the center of the pore. The four pores show the same 
feature, but the narrower channels restrict fluid back-flow. 
The consequence is that pressure builds up in the hot bulk 
regions, which is the thermo-osmotic effect. Preliminary data 
for AP and Dp for this one-component fluid are given in 
Table 2. These results are for a one-component fluid, so there 
is no Soret effect here. We will present results for various 
porosities and permeabilities for one- and two-component 
fluids at IMT15. 


System layout 


Velocity fields 


@ a 


40 


z* 20 20 
10 10 
C) 0) 
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Figure 1. Top: Side view of two slit pore configurations used in this 
work, one wide pore (left) and four narrow pores (right). The blue 
regions are pore walls of Lennard-Jones/spline particles in a solid 
state. The red regions contain the same type of particles in a one- 
component liquid state. A temperature gradient is set up in the fluid 
between the hot and cold bulk compartments. Middle: Convective 
flow fields in the pores as illustrated by the velocity vectors. The 
velocity is towards the hot bulk along the walls and in the opposite 
direction in the pore centers. Bottom: x-components of the fluid 
velocity in the left pores as function of z. 
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Table 2. Preliminary results for AP* and D5 (in Lennard-Jones 
units) for the two pores shown in Figure 1. The thermo-osmotic 
effect in the wide pore is very small because of the low resistance 
to fluid backflow in the center of the pore. The pressure build-up on 
the hot side is significant in the system with narrow pores. 


Case AP* (hot-cold) Dp 
0.016 
0.174 


Wide pore 


Narrow pores 
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Introduction 


Motion of solutes in a solvent induced by a thermal gradient 
is termed as thermophoresis. Ratio of the established 
concentration gradient to the temperature gradient is 
quantified using Soret coefficient St (Kéhler et al., 2016). 
Since this phenomenon is very sensitive to the nature of 
solute-solvent interactions, it is used asa tool for quantifying 
biomolecular interactions, especially _ protein-ligand 
interactions (Niether et al. 2020). Although, a change in the 
thermophoretic behaviour of protein once the ligand binds is 
often attributed to the change in the hydration layer, 
underlying microscopic physical effect is not understood. To 
gain deeper insight into the interactions involved, we 
investigate whether the non-equilibrium coefficient measured 
can be related to equilibrium properties. The same, 
thermophoretic data measured using thermal diffusion forced 
Rayleigh scattermg (TDFRS) (which is a non-equilibrium 
process) are compared with the thermodynamic data obtained 
by isothermal titration calorimetry (ITC) (which is an 
equilibrium method). To connect these parameters, we start 
from an early work by Eastman (Eastman 1928) which 
connects Srto Gibb’s free energy AG as follows 
S,= (1) 
kpT aT 

Integrating Eq. | with temperature gives us access to a 
relation between Sr and AG for the individual compounds of 
the system (free protein, free ligand, and complex). These 
individual contributions calculated are then used to establish 
arelation between ITC and TDFRS measurements, as shown 
in Fig. 1. We hypothesize that AGr nigh can be calculated 


from the free energy change at a low temperature AG, 
measured by ITC and the differences in AAG corresponding 
to two temperatures for the individual components probed by 
TDFRS using the following equation 


AGrnign = AG rio + SAGaz — AAG, — AAG, — (2) 


Figure 1: Schematic illustration of the calculation of AG and AAG 
from ITC and TDFRS, respectively. 


Choice of systems 

To test our hypothesis, EDTA-CaCh reaction in MES buffer 
is used as reference system. This is a well-known chelation 
reaction which is used as a validation standard for ITC 
measurements (Rafols et al. 2016). Once the validation is 
done, Eq. 2 is tested for protein Bovine Carbonic Anhydrase 
I (BCA I) with two different ligands which belong to the 
group of arylsulfonamides. In our study, we used 4 
fluorobenzenesulfonamide (4FBS) and 
pentafluorobenzenesulfonamide (PFBS). We carry out 
TDFRS and ITC measurements for these systems overa wide 
temperature range between 20°C to 45°C. 


Results 


1. TDFRS measurements 
1.1 EDTA-CaCh 


St of EDTA shows an abrupt drop between 25°C and 30°C. 
Although, there are several literature reports, which report a 
change in properties of several systems in the presence of 
EDTA in a similar temperature range, there exists no clear 
explanation for this behavior (Sugiyama ef al. 2014, 
Minkevich et al. 2006, Yilmaz et al. 2011). Also, the 
temperature sensitive polymer poly(N-isoproplacrylamide) 
(PNiPAM) in water (Kita et al. 2005), which is influenced by 
the hydrogen bonds, shows an abrupt variation in Sr with 
temperature. Distinctive thermodiffusive properties exhibited 
by EDTA and the complex might have the same origin as in 
case of PNiPAM as it happens in similar temperature range. 
This also supports the hypothesis that a change in hydration 
hasan effect on thermophoretic behavior. 


1.2 Protein-ligand system 


- ow" 
@- & @ Pres 
k= ee é @ 4FBS 
L Lee @ BCAI4FBS 4 
é-- @ BCAI-PFBS 
@ BCA! (pH7.4) 
X Buffer 


Oolx % “x ™& X*& ™& | 
20.25 30 35 40 45 
temperature / °C 


Figure 2: Temperature dependence of St for BCA I, 4FBS, PFBS, 
buffer and complexes formed. 


Temperature dependence of St for protein-ligand systems is 
shown in Fig. 2. From TDFRS measurements, it can be 
concluded that the hydration shells of both complexes formed 
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are very similar, but different from free BCA I. For both 
protein—ligand systems the difference between St of free 
protein and complex decreases with temperature, so that it is 
almost negligible at high temperatures. This is an indication 
that the binding of both ligands should become weaker with 
increasing temperature. A larger change in Sr with 
temperature for free BCA I indicate a higher hydrophilicity 
for the free protein in comparison to the protein-ligand 
complex formed. It has been reported that arylsulfonamides 
bind to the active site of different variants of carbonic 
anhydrase protein by replacing a water molecule from the 
active site (Krishnamurthy et al. 2007, Olander et al. 1973, 
Abbate et al. 2002). This implies that the complex is less 
hydrophilic than the free protein which is confirmed by 
TDFRS measurements. 


2. Validation of relation between St and AG 


The calculated values listed in Table 1 correspond to Thigh = 
30°C and Tiow = 20°C. Mathematical expression is first 
applied to the EDTA-CaCh system. With the St values of 
EDTA, CaCl2 and the complex measured at Thighand Tiow, we 
have access to the change in Gibb’s energy (AAG) of the 
individual components. Later, the same procedure is used for 
protein-ligand systems. 


Table 1: Comparison of AG that has been calculated and that has 
been measured with ITC 
AGealculated AGitc 


System 
(kJ/mol) (kJ/mol) 
EDTA-CaChk -36.5+1.2 -36.4+0.8 


-38.2+1.5 


BCA I-4FBS -40.541.1 -40.4+41.3 


BCA I-PFBS -39,943.9 


For all the three systems that have been studied, AG that has 
been calculated agrees with the measured values within the 
error limits. 


Conclusions 


One of the main goals of this work was to investigate whether 
it is possible to connect thermodynamic parameter obtained 
with ITC to thermodiffusion parameter measured with 
TDFRS. A mathematicalrelation connecting Stand AG was 
derived which is tested for three different systems; low 
molecular weight reference system, EDTA—CaCh and the 
protein BCA I with two ligands 4FBS and PFBS. For all the 
temperature pairs that have been studied for aforementioned 
systems, AG values calculated agreed with the measured 
values within the error limits. This implies that St measured 
at different temperatures can be used to predict the AG value 
at particular temperature. This newly developed connection 
can be utilized to open promising gates in the accurate 
acquisition of the underlying binding and molecular 
dissociation mechanisms from the studied systems. 
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Introduction 


Phase Change Materials (PCMs) are of interest in thermal 
management for their ability to store (and release) large 
amounts of energy, in the form of latent heat, near the phase 
change temperature. Space exploration is an especially fitting 
area for PCM applications because significant variations in 
temperature can result from cycles of radiation exposure or 
the heat generated by onboard systems. Although organic 
PCMs are a stable and attractive choice for many situations, 
their effectiveness can be limited by low thermal 
conductivity. 


One straightforward strategy for improving performance is to 
modify the (effective) thermal conductivity of the PCM 
device by adding a more conductive material — metallic 
structures or nanoparticles, for example. Alternatively, 
encouraging convective flows in the liquid phase can 
significantly enhance the heat transfer rate. On ground, 
gravity can promote natural convection if the PCM device has 
a suitable orientation. In microgravity environments, the 
thermocapillary effect can similarly drive convective flow 
and has been proposed as an alternative strategy since, in the 
presence of a free surface, the temperature gradient inherent 
to PCM operation will naturally generate surface tension 
gradients. 


In recent years, the topic of thermocapillary-driven melting of 
PCMs in microgravity has been the subject of growing 
research. The seminal parabolic flight experiments of 
Ezquerro et al. (2020) were the first to confirm the potential 
of the thermocapillary effect to enhance heat transport during 
melting in microgravity. Since then, severalnumerical studies 
have looked at this problem from the perspectives of heat 
transport (see, e.g., Salgado Sanchez et al. 2020, Varas et al. 
2021, and reference therein) and pattern selection during the 
phase change. In fact, the relevance of this research line for 
space exploration was recently recognized by the European 
Space Agency with the approvalof the ‘Effect of Marangoni 
Convection on Heat Transfer in Phase Change Materials’’ 
(MarPCM) experiment (Porter et al. 2023) for its future 
execution on board the International Space Station. 


In this work, we review the state of the art available in the 
literature and discuss different strategies that can be 
implemented in microgravity to improve the performance of 
thermocapillary-enhanced PCM devices: selecting an 
adequate geometry, using (nano-enhanced) NePCM, or a 
combination of both. 


Selecting PCM geometry 


The types of PCM containers that been analyzed in 
thermocapillary-driven melting can (to the best of our 
knowledge) be placed into three general categories: 
quadrilateral (rectangular and trapezoidal cross-section), 


cylindrical and spherical (i.e., liquid bridges and drops when 
melted). 


In rectangular containers, the effect of thermocapillary flow 
was shown to depend crucially on aspect ratio [T = L/H and 
the Marangoni number 
ae yL =m 
La 

which quantifies the relative importance of thermocapillary- 
and conductive transport during the phase change. Here, y is 
the thermocapillary coefficient, L is the containerlength, H 
is itsheight, AT is the applied temperature difference, and yu 
and a@ are the viscosity and thermal diffusivity of the liquid 
phase, respectively. The enhancement factor G = Tye¢/T ya; 
defined as the ratio between the melting time with purely 
conductive transport T,¢¢ and with the thermocapillary effect 
Tma> fanges from approximately 2 to 25 (Salgado Sanchez et 
al. 2020). For fixed AT, there exists an optimal IT" that 
maximizes this enhancement. 


With rectangular containers, the phase change often includes 
a final stage that is dominated by conduction, where the PCM 
near the cold boundary melts slowly. Recently, the use of 
trapezoidal containers was proposed to accelerate the melting 
process. By inclining the cold lateral wall, in particular, 
melting can be substantially speed up with respect to the 
rectangular case, with an optimal reduction of Ty, by a 
factor of 3 in the limiting case of a right triangle (Borshchak 
Kachalov et al. 2022). 


In liquid bridge geometry, on the other hand, the contribution 
of the thermocapillary effect during the phase change was 
shown to be approximately 50% greater than in the 
rectangular case (Varas et al. 2021). Note that, for a fixed heat 
storage capacity (i.e., PCM volume), this geometry increases 
the ratio of free surface area to PCM volume compared to the 
most trapezoidal (orrectangular) containers. With this criteria, 
the optimal geometries are spherical shapes, for which the 
free surface driving the thermocapillary flow is maximized 
for a given PCM volume. 


Using NePCM with thermocapillary effects 


The use of dispersed nanoparticles has also been proposed for 
improving PCM performance with the (so-called) nano- 
enhanced phase change materials (NePCM) (Madruga et al. 
2021). 


The addition of nanopartciles modifies the thermo-physical 
properties of the pure PCM, increasing thermal conductivity 
(diffusivity) and the viscosity of the liquid phase, and slightly 
reducing the heat storage capacity; note that part of the 
original PCM volume is taken by particles that do not undergo 
a phase change. The extent of this modification depends on 
the concentration C and size of the particles, and their 
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thermo-physical properties. Compared to the response of the 
pure PCM at a given AT, the associated NePCM will 
experience a smaller Marangoni number and possess lower 
latent heat. 


However, the interplay between these two opposite effects in 
rectangular geometry could result in an overall enhancement 
ranging from 1 to 10%, depending essentially on 
nanoparticles properties and concentration C, and on I and 
AT. 


Combining strategies 


For applications, one may consider combining the 
aforementioned strategies: selecting an adequate (optimal) 
geometry for an NePCM. In this case, the contribution ofeach 
strategy would be complementary. A rectangular container 
with NePCM can provide values of G ranging from 1% to 
10% greater than with the pure PCM. 


In cylindrical geometry, the factor of 1.5 achieved by the 
greater free-surface-area-to-volume ratio would be increased 
by an additional factor of approximately 1.1, resulting in a 
final enhancement of roughly 1.65. Compared to the case of 
purely conductive transport, the overall enhancement would 
include contributions from: the thermocapillary flow, with a 
G factor on the order of 10, the container geometry and the 
NePCM, which provide additional improvements of 50% and 
10%, respectively. Overall, an NePCM in a liquid bridge 
geometry would melt roughly 17 times faster than the same 
configuration of a pure PCM without a free surface. 


Conclusions 


The present and future of space exploration requires improved 
thermal control systems. Organic phase change materials 
(PCMs), with their large latent heat, suitable melting 
temperatures and chemical stability, represent an attractive 
option for passive thermal control. 


Recently, the use of thermocapillary flows has been shown to 
be a viable strategy for improving the low heat transfer rate 
of organic PCMs, whose performance is otherwise hampered 
by long melting/solidification cycles. In this work, we 
summarize the basic options that have been proposed to 
accelerate the thermocapillary-driven melting of PCMs in 
microgravity: an appropriate selection of container geometry 
and the use of dispersed nanoparticles to increase the effective 
conductivity of the PCM. Compared to the case of pure 
conduction, the complementary improvments from the 
thermocapillary effect, container geometry and the addition 
of nanoparticles can increase the melting rate by an order of 
magnitude. 


Tarragona (Spain), May 2023 
Acknowledgements 


This work was supported by the Ministerio de Ciencia e 
Innovacion under Project No. PID2020-115086GB-C31, and 
by the Spanish User Support and Operations Centre (E- 
USOC), Center for Computational Simulation (CCS). 


References 


J. M. Ezquerro, P. Salgado Sanchez, A. Bello, J. Rodriguez, V. 
Lapuerta and A. Laverén-Simavilla, Experimental evidence 
of thermocapillarity in phase change materials in 
microgravity: measuring the effect of Marangoni conbvection 
in solid/liquid phase transitions, International 
Communications in Heat and Mass Transfer, 113, 104529 
(2020) 


P. Salgado Sanchez, J. M. Ezquerro, J. Fernandez and J. 
Rodriguez, Thermocapillary effects during the melting of 
phase change materials in microgravity: heat transport 
enhancement, International Journal of Heat and Mass 
Transfer, 163, 120478 (2020) 


R. Varas, P. Slagado Sanchez, J. Porter, J. M. Ezquerro and V. 
Lapuerta, Thermocapillary effects during the melting in 
microgravity of phase change materials with a liquid bridge 
geometry, International Journal of Heat and Mass Transfer, 
178, 121586 (2021) 


J. Porter, A. Laveron-Simavilla, M. M. Bou-Ali, X. Ruiz, F. 
Gavalda, J. M. Ezquerro, P. Salgado Sanchez, U. Martinez, D. 
Gligor, I. Tinao, J. Gomez, J. Fernandez, J. Rodriguez, A. 
Borschak Kachalov, V. Lapuerta, B. Seta, J. Massons, D. 
Dubert, A. Sanjuan, V. Shevtsova and L. Garcia-Fernandez, 
The “Effect of Marangoni Convection on Heat Transfer in 
Phase Change Materials” experiment, submitted to Acta 
Astronuatica, under review. 


A. Borshchak Kachalov, P. Salgado Sanchez, U. Martinez, J. 
Fernandez and J. M. Ezquerro, Optimization of 
thermocapillary-driven melting in trapezoidal and triangular 
geometry in microgravity, International Journal of Heat and 
Mass Transfer, 185, 122427 (2022) 


S. Madrugaand C. Mendoza, Heat transfer performance and 
thermal energy storage in nano-enhanced phase change 
materials driven by _ thermocapillarity, International 
Communications in Heat and Mass Transfer, 129, 105672 
(2021) 


66 


15" International Meeting on Thermodiffusion (IMT15) 


Tarragona (Spain), May 2023 


On the potential of thermodiffusion as means for large-scale desalination 


Shuqi Xu! and Juan F. Torres'” 


' School of Engineering, Australian National University, Canberra, ACT 2601, Australia 
* Email: felipe.torres@anu.edu.au 


Introduction 

Thermodiffusion (TD) has found use in few yet essential 
engineering applications, such as enrichment of nuclear 
material (Nier, 1989) and analysing protein binding (Wienken 
et al., 2010). Despite TD being a phenomenon first discovered 
by Soret in 1879 through a desalination process (Soret, 1879), 
it has never been used as means of desalination because it is 
a weak phenomenon that may not meet the requirements of 
purity (<500 ppm) and yield (+100 m?/day) for large-scale 
desalination. A facile method for removing salt ions from 
water could be key in alleviating water scarcity, which is one 
of the most pressing issues of our time (World Bank, 2019). 
Conventional desalination methods such as distillation and 
reverse osmosis are either very energy intensive or have 
membrane fouling issues. To the best of our knowledge, there 
are only two theoretical papers exploring the possibility of 
thermodiffusive desalination (TDD) including our 
preliminary work (Abdeljabar, 2011; Xu et al., 2022). TDD 
is appealing because it is a membrane-free method that can be 
driven by low-grade thermal energy from industrial waste 
heat, solar energy or the surrounding environment. Here, we 
explore the possibility of scaling up TDD and reflect on its 
challenges and prospects. First, we provide experimental 
results showing that TD is able to achieve a concentration 
drop (Carp) of 450 ppm and 700ppm from a NaCl 
concentration of 30,000 ppm (seawater) and 60,000 ppm 
(brine), respectively, when passing the solution through a 
channel. Since Cuarop is relatively small after a single pass, 
multiple passes are implemented to increase it. However, the 
recovery rate for that concentration drop is reduced to 6% (i.e. 
the percentage of low salinity water extracted from a high 
salinity water volume). A cascaded microfluidic structure that 
can potentially overcome this challenge and scale up TDD is 
proposed and theoretically assessed. 


Single-pass TDD 

The concept of TDD is depicted in Figure 1a. A laminar flow 
of an initially homogeneous saline solution passes through a 
rectangular channel with its upper and lower walls kept at 
constant high and low temperatures, respectively. NaCl is 
thermophobic in water above 12°C (Caldwell, 1973), thus the 
salt migrates to the cold side in our thermodiffusion 
experiment. Together with the thermal stratification, solute 
density stratification stabilises the flow. TDD is potentially a 
low-energy consumption desalination methods as it does not 
involve phase change (as in distillation) and can be driven by 
low-temperature heat from the sun or industrial processes. 
Phase-shifting interferometry (PSI) (Torres eft al., 2013) is 
used to accurately measure the concentration difference AC 
between the top and bottom streams, which come from an 
equally bifurcated flow at the end of the rectangular channel. 
As shown in Figure 1b, we control the flow rate while 
degassing the saline solution to avoid bubble-induced mixing. 
A smaller channel height means faster thermodiffusive 


separation but also larger heat flux. Discussion on the channel 
dimensions and other control parameters can be found in our 
previous numerical work (Xu et al., 2022). 


AC between the two streams was precisely measured using 
PSI and the results are shown in Figure 2. ATmeas is obtained 
from the mean measured temperatures from thermocouples 
placed along the channel, while A7set refers to that set in the 
two water baths. The effect of A7imeas on separation is shown 
together with CFD results assuming a fully-developed 
velocity field. For the same mean temperature Tinean, we could 
confirm that the larger the AZmeas, the larger the AC. In 
addition, since the Soret coefficient Sr is larger at higher 
temperatures (Caldwell, 1973), AC increases with Tinean for 
the same ATmeas. The effect of flow rate was also explored. 
For a flow rate less than 5 mL/min, saline water stays in the 
1 mm-high channel for more than 2 min, achieving a near 
maximum AC. Larger flow rates drop AC as there is not 
enough residence time of the solution in the channel for 
thermodiffusion to fully separate the NaCl species before 
reaching the bifurcation. 


(a) 


PSI 


Low 


<a concentration 


\gootine 


®\ High 
concentration \ 


Figure 1: Concept demonstration and experimental setup of the 
TDD unit. (a) Seawater flows through the TDD channel with a 
positive temperature gradient across its height, then separates into 
two streams bifurcated at the exit. AC between bifurcated water 
streams is indicated by digital interferometry. Here, the top-stream 
sample (red square) was injected into a test cell half-filled with 
bottom-stream sample (blue square). The fringes in the green 
square indicate AC when unwrapped. (b) The flow rate of the saline 
water in the channel is controlled by a peristaltic pump. Saline 
water is degassed before entering the channel. The fluid path for the 
saline water is indicated by orange lines. The channel is 500 mm 
long, 20 mm wide, and | mm high. The collection bottles for the 
two streams are at the same height, ensuring equal flow rate. Water 
circulating from two separate water baths (blue and red arrows) 
creates the AT. Six pairs of thermocouples were distributed 100 mm 
apart on both sides of the channel. 
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Figure 2: Single-pass TDD quantified by the AC between the two 
water streams with different Tinean and ATineas. 


Multiple passes 

As shown in the single-pass case (Figure 2), the maximum 
AC between the top and lower flow streams is ca. 900 ppm 
for AT set = 60 K (ATimeas =37 K). Above ATset = 60 K, the heat 
exchange between the two water baths was too large, resulting 
in a marginal increase in A7meas. Note that 900 ppm is 3% of 
the initial seawater concentration (relative value). Moreover, 
the concentration drop Carop of the top stream relative to the 
inlet concentration is only half of that value (i.e. 450 ppm). 
This small drop in concentration is not enough to obtain 
freshwater water useful for agriculture (accounting for ca. 
70% of water use worldwide), which requires 95% relative 
concentration drop. To amplify the thermodiffusive 
separation, a simple way is to accumulate the top low- 
concentration stream in a container and then pass it again 
through the channel. This experiment was done for both 
seawater and brine, with results shown in Figure 3. Carop 
increases with the number of passes, demonstrating that TDD 
can be enhanced. Furthermore, it is noticed that Carop is larger 
in brine than in seawater. This is because thermodiffusion 
increases with concentration (reaching a maximum at C = 0.5 
as per thermodiffusion mass flux equation (Xu et al., 2022)). 
This result suggests that TDD may be a potential pre- 
treatment method when using brine as the RO water feed 
because RO is more energy intensive when C is high. The 
significant drawback of this multiple-passes TDD is that the 
volume of the desalinated stream drops exponentially. In the 
last pass, the volume (i.e. recovery rate) was too small that the 
contamination in the channel and the peripherals diminished 
the concentration drop due to TDD. 


This multiple-passes TDD is not efficient because the ‘high’ 
concentration solution from a pass may have a lower 
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Figure 3: Multiple-pass TDD of NaCl aqueous solution. The 
concentration drop is plotted as a function of number of passes 
recirculated for two different initial concentrations (30,000 
and 60,000 ppm). 
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concentration than the initial saline water (seawater or brine) 
in the system, but it is nonetheless discarded. A scalable 
microfluidic thermodiffusive separation device known as the 
Burgers cascade was proposed for enhancing thermodiffusive 
separation in gases (Saiki et al., 2020). It is a two-dimensional 
device made up of many small cells, each similar to a mini 
TDD channel. Different species in seawater have different Sr 
(Caldwell and Eide, 1981). If assuming seawater behaves as 
a binary mixture and has a bulk Str of 0.005 K", the single- 
pass TDD CFD model can be applied to this device. 
Concentration in each cell was theoretically modelled and the 
low concentration of 1000 ppm was predicted for a recovery 
rate of 10%. Future work will experimentally validate the 
assumption of bulk thermodiffusion in multicomponent 
electrolytes and improve the Burgers cascade design to 
maximise the recovery rate and concentration drop. 


Conclusions 

It was experimentally shown that thermodiffusion can drop 
the concentration in a binary approximation of seawater (or 
brine) by around 450 ppm (or 700 ppm) in a single pass, 
which is only about 3% of the initial concentration. However, 
at least 95% of salt needs to be removed to obtain fresh water. 
In an attempt to amplify Carp, multiple-pass TDD was used 
and Carop > 1000 ppm was achieved, but an exponential drop 
in recovery rate would render TDD unfeasible. To address this 
problem, a scalable microfluidics device for TDD is proposed 
an analysed theoretically. 
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Introduction 


Fluid systems containing bubbles are widely observed in 
natural phenomena such as breaking waves (Deike et al. 
2016} and a wide range of engineering applications including 
bubble column reactors (Kantarci et al. 2005) and heat 
exchangers (Sadighi 2014). Momentum and_heat/mass 
transfer between the bubbles and the surrounding fluids are 
common phenomena in many industrial processes including 
chemical reactors based on bubbly flow (Schluter et al. 2021), 
multiphase cooling systems (Suranjan et al. 2009), 
wastewater treatment (Shirai et al. 1999), metal refining 
(Wang et al. 1996) and also in nature such as bubbly plumes 
under the oceans (Boufadel et al. 2020). 


Generally, buoyancy-driven bubbly flows are utilized to 
improve heat/mass transfer behavior in industrial 
applications. Heat/mass transfer across the bubble's interfaces 
can be affected by the flow behavior inside and around the 
bubbles, flow characteristics, topological changes in the 
bubbles, distributions of the bubbles, interactions between the 
bubbles, and bubble coalescence or breakup. To measure the 
heat/mass transfer rate, advection and diffusion within each 
phase and also transfer across the interface should be 
considered. Prediction of the heat/mass transfer can be a 
challenging problem due to the bubbles interactions, 
deformable interface, and differences between the diffusion 
coefficients in each phase (which could be of orders of 
magnitude). Interactions between the bubbles is one of the 
most important issues in bubbly flows which can lead to 
important physical processes occurring in the system such as 
coalescence and breakup. For example, coalescence can 
affect the performance of the system by changing the 
heat/mass transfer rate caused by altering the interfacial area. 
The evolution, flow dynamics and scalar transport in a system 
comprised of only two bubbles can be of help to better 
understand the nature of bubbly flows and improve the design 
of gas-liquid equipment. In this context, it is necessary to have 
sufficient knowledge of the bubbles dynamics, especially 
interactions between the bubbles, to predict the heat/mass 
transfer rate in bubbly flows. 


A numerical simulation was performed to study the motion 
and interaction of two in-line bubbles rising in Newtonian and 
shear-thinning fluids and mass/heat transfer across the 
bubbles interface was investigated. A three-dimensional 
Volume-of-Fluid method (VOF) was used to capture the 
interface and the Carreau model was used to describe the 
rheological behavior of the shear-thinning fluid. The 
dynamical behavior of a single air bubble rising in water 
depends on the Galilei (Ga) and Bond (Bo) numbers shown 
in a Ga-Bo phase plot in Figure 1. As can be seen, there are 
five different regimes including I. axisymmetric, II. skirted, 


IU. oscillatory, IV. peripheral breakup, and V. central 
breakup. Here, the dynamical behavior of a pair of bubbles 
and mass transfer across the bubbles interface in each regime 
have been studied both in a Newtonian fluid and a shear- 
thinning fluid. Moreover, the effects of the rheological 
properties such as flow index (n) and inelastic time constant 
(A) on the mass transfer across the bubbles interface and on 
the bubbles dynamics have been investigated for a fixed set 
of (Ga,Bo). 


Figure 1: Ga-Bo phase plot of an air bubble rising in water. 


Conclusions 


It is found that the dynamical behavior of a pair of bubbles 
can be different from the one of a single bubble and 
decreasing the flow index leads to increasing local Ga. \ has 
a small effect on the mass transfer rate while the oscillatory 
behavior of the bubbles intensifies by increasing i. 
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Introduction 

Non-equilibrium fluctuations (NEFs) are a model for 
investigating the diffusion processes in complex fluids. 
Shadowgraphy is a robust near-field scattering technique with 
a simple configuration and is highly sensitive to NEFs. 
Besides allowing the measurement of the transport properties 
in complex fluids (Croccolo et al. 2006a, Croccolo et al. 
2012), shadowgraphy allows probing the effect of diffusion, 
gravity (Croccolo et al. 2006b) and confinement (Giraudet et 
al. 2015) on the dynamics of NEFS and the presence of 
propagating waves due to the coupling of the different 
relaxation modes (Croccolo et al. 2019; Garcia-Fernandez et 
al. 2019). The need for developing the two-wavelength 
shadowgraphy is dictated by the impossibility of separating 
the two solutal relaxation modes in a ternary molecular 
mixture with a classical one-wavelength shadowgraph 
(Bataller et al. 2016). Indeed, in such mixture, the eigenvalues 
of the diffusion matrix are almost identical, implying that the 
relaxation times of the NEFs of the two independent 
concentration fields cannot be distinguished. In such case, it 
is impossible to differentiate the two eigenvalues of the 
diffusion matrix. Using the same approach as the one 
implemented in the two-wavelength selectable optical digital 
interferometry (SODJ) of the DCMIX project (Mialdun et al. 
2020) to separate the two independent concentration fields of 
a ternary mixture, a two-wavelength shadowgraph would 
allow, in principle, to separate the contributions of the two 
concentration fields to the shadowgraph signal. Indeed, the 
same approach can be used to separate NEFs of temperature 
and concentration fields in a binary mixture subjected to a 
temperature gradient. The binary system is a simple and 
effective test for the two-wavelength shadowgraphy 
technique. In this case, the relaxation modes of temperature 
and concentration are well distinct in time and can therefore 
be also separated with a classical shadowgraph by fitting the 
structure function of the NEFs with two exponential decays. 
In this work, we present the principles of two-wavelength 
shadowgraphy and demonstrate the technique working 
principle by performing a thermodiffusion experiment on an 
equimassic binary mixture of toluene/n-hexane subjected to a 
temperature difference of 20 K. 


Working principle of two-wavelength shadowgraphy 

The shadowgraph consists of two superluminescent diodes (A 
= 675+ 13 nm or A= 405+ 13 nm) connected to a single-mode 
optical fiber that delivers a weakly coherent beam collimated 
by an achromatic lens. The intensity of the light scattered 
strikes the CCD sensor of an s-CMOS camera placed in the 
near field of the thermodiffusion cell containing the sample. 
Owing to their random nature, NEFs are characterized by an 
autocorrelation function and static structure factor. These 
quantities can be calculated theoretically by fluctuating 


hydrodynamics (Ortiz De Zarate et al. 2006) and are directly 


related to the static structure factor and autocorrelation 
function of the refractive index fluctuations measured by 
shadowgraphy. The case of an ordinary binary mixture (with 
a large Lewis number) subjected to a temperature gradient in 
the Boussineq approximation (pressure fluctuations are 
neglected) is considered here. With a large Lewis number, all 
couplings between the temperature and concentration 
fluctuations are neglected, at least for large enough wave 
numbers, and the autocorrelation function of the refractive 
index fluctuations for a given wavelength A can be written 
as: 
2 2 

5x(q, At) = (F) Seclq, At) + (F)_ Sera. At) (1) 
where q is the wave number, i.e. the modulus of the wave 
vector. The derivatives of the refractive index with respect to 
the temperature and concentration depend on the wavelength 
of the probe beam. Using shadowgraphy, the wave vector at 
which the structure factor or autocorrelation function is 
measured does not depend on the wavelength of the probe 
beam or scattering angle, as is the case in far-field scattering. 
Therefore, for the same wave number gq, the autocorrelation 
function of NEFs with probe beams of different wavelengths 
can be easily measured. At large wave numbers, any coupling 
between different modes, which is responsible for the 
appearance of the propagating modes, disappears. The 
concentration and temperature signals can be easily separated 
if the optical contrast factor matrix is known: 


Scc(q, At) (), (7). 5,,(q, At) 
7 : (2) 
ey i) Vso 


Result and discussion 

From the shadowgraph images, one can thus calculate the 
structure function (figure 1. a), and the static structure factor 
(figure 1. b), from which the autocorrelation functions (figure 
l.c) of the NEFs at the two wavelengths are calculated. 
Having previously measured the optical contrast factors by 
means of a refractometer, the autocorrelation functions of the 
temperature and concentration NEFs can be separated (figure 
1. d). On the one hand, the individual analysis of the structure 
functions makes it possible to determine the relaxation times 
of the NEFs of temperature and concentration, as well as their 
intensities. On the other hand, the analysis of the 
autocorrelation functions of the NEFs of concentration and 
temperature, separated optically, makes it possible to 
determine the relaxation times distribution of the fluctuations 
as well as their intensities. Preliminary analysis shows that the 
results obtained by the two procedures correspond for the 
majority of the wavenumbers. Further analysis is in progress, 
and the results will be presented at the IMT15. 


Spr (q, At) 
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Figure 1 (a) Structure functions for both wavelengths and for g = 
118 cm. (b) Static structur factor for both wavelengths. The vertical 
dashed line stands for g = 118 cm (c) Autocorrelation functions 
extracted from (a) and (b) and for g = 118 cm. (d) Extraxted 
autocorrelation functions of the temperature and concentration NEFs 
as obtained by applying Eq. 2 to the data of panel (c). The 
distributions of the relaxation times of the temperature and 
concentration fluctuations are determined from the separate signals. 
Analysis of the distribution of the relaxation times makes it possible 
to determine the thermal and mass diffusion coefficients. 


The procedure for separating the temperature and 
concentration modes in a binary system is robust enough 
despite a large value of the condition number of the contrast 
factors matrix (>400). While the analysis involved in SODI is 
based on fitting through a series of exponential decays, the 
Fourier-based analysis of Shadowgraphy is based on simple 
exponential decays. This actually increases the applicability 
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of the Shadowgraph technique to a larger range of condition 
numbers. The large condition number obtained here, is due to 
the fact that the thermal contrast factor varies very little with 
the wavelength. This result validates the technique and 
promises even better results when separating solutal modes in 
a ternary mixture. 


Conclusions 

In this work, we make a first demonstration of the faiseability 
of the two-wavelength shadowgraphy technique to separate 
the signals of uncoupled fluctuating fields in a binary or 
ternary liquid mixture. First, we tested the ability of the 
technique in separating the temperature and concentration 
fields in a binary mixture subjected to a temperature gradient. 
The results obtained make it possible to validate the technique 
and demonstrate its robustness. The next step of our work will 
be to test the technique on a ternary polymer mixture and then 
on a ternary molecular mixture. 
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Introduction 

Fluctuations occur spontaneously in a fluid and are amplified 
by a gradient of at least one thermodynamic variable 
(concentration, temperature or pressure). The resulting non- 
equilibrium fluctuations (NEFs) can be studied by means of 
several optical techniques, such as Dynamic Shadowgraphy, 
Dynamic Schlieren and Differential Dynamic Microscopy 
(DDM). The spatio-temporal evolution of NEFs, can be 
represented by the structure function (Croccolo et al. 2007, 
Croccolo et al. 2012, Schulz-DuBois and Rehberg 1981) over 
a wide range of wave numbers, as calculated by means of the 
differential dynamic algorithm (DDA). The structure function 
is obtained using a series of images obtained by one of the 
optical techniques cited above (Croccolo et al. 2007, 
Norouzisadeh et al. 2021). At each wave number, the structure 
function can be analyzed after providing a theoretical model 
(fitting equation) to extract the transport properties of the 
fluid under study, such us the mass diffusion coefficient, 
critical cut-off vector, and eventually from those values, the 
Soret and thermodiffusion coefficients can be derived. The 
applicable model varies depending on the nature and 
conditions of the experiment as well as the wavenumber. The 
best model should be used for each wave number in order to 
obtain a suitable fitting, which can be a challenge for the 
analysis of large amounts of data, especially when different 
physical phenomena dominate at different wavenumbers. 

In this study, we take advantage of artificial intelligence (AI) 
to analyse the experimental structure function. We developed 
a program based on neural networks that chooses the best 
model to fit the structure function for each wave number thus 
allowing an automatic, fast and efficient analysis independent 
of the amount of data to be treated. 


Structure function analysis 

Fitting theoretical models to experimental data is important in 
all scientific fields. It aims at estimating unknown parameters 
of models from data using a variety of statistical 
methodologies. In our case, the experimental structure 
function should be fitted with a theoretical model to extract 
the transport properties of the fluid under study. The basic 
theoretical equation of the structure function is given by: 

SF = 2{T(q)S(q)[1 — ISF(q, At)] + BG(q, At)}, (1) 
where q and At represent the wave number and the 
correlation time, respectively. T(q) is the optical transfer 
function of the shadowgraph, S(q) the static power 
spectrum of the fluid density fluctuations and ISF(q, At) is 
the Intermediate Scattering Function. The latter corresponds 
to the dynamic part of the structure function that can be 
modelled by a sum of exponential decays: 

At 


ISF(q, At) = Yin, aexp(— 7); (2) 
where n is the number of relaxation modes in the system (e.g. 
concentration and temperature), a; are the relative 


amplitudes of the different modes such as }};a; = 1 and the 
T;(q) represent the decay times. 

At small wave number, in the presence of propagating modes, 
a sinusoidal term needs to be included to describe the 
oscillations in the SF originated by the presence of 
propagating modes (Croccolo et al. 2019): 


_ = At az 
ISF (q, At) = a, exp( a + coslot@l x cos[OQ(q)At + 
t 
o(q)]exp(—p). (3) 


where 1 and @ depend on q and represent the oscillation 
frequency and phase. 
In Eq. 1, BG(q, At) represents the optical background that 
can take different forms depending on the experimental 
conditions. Here, we take two possible forms, constant or 
quadratic: 
BG,(q, At) = A(q), (4) 

BG,(q, At) = A(q) + B(q)At?, (5) 
where A describes the experimental noise, while B is related 
to the stability of the background itself. 


Automatic analysis 

The experimental structure function can be fitted using 
different models, depending on the investigated system. 
Manual or non-automatic analysis of a large amount of data 
is time consuming and may lack accuracy. Indeed, in a 
thermo-diffusion experiment, the propagation modes are 
expected at small wave numbers and both the concentration 
and temperature modes are expected at different wave number 
ranges. In addition, the optical background can take different 
shapes from one wave number to another. In this study, we 
show how we can leverage the power of the deep learning, 
especially convolutional neural networks (CNN) (Albawi et 
al. 2017) to automatically analyze the experimental structure 
functions. Six possible fitting equations were 
considered, JSF in the propagation modes, JSF with one or 
two exponentials (n=1 or n=2), and constant or 
quadratic background. The output coefficients are: T(q)S(q), 
ai(q), t1(q), A(q) and B(q). 

As shown in Fig. 1, our approach consists in passing the 
image pixel values of the structure function shape through 
several convolutional layers looking for patterns and creating 
feature maps that represent the inputs of a fully connected 
artificial neural network (ANN), to finally obtain the best 
theoretical model for analyzing the structure function. To do 
so, three steps were followed. First, a “learning phase” allows 
creating learning classes containing different representative 
structure function shapes. Homogeneity within each class 
should be maximized while homogeneity between classes 
should be minimized in order to reach high accuracy results. 
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Figure 1: A typical CNN architecture scheme applied on recognizing the more adequate fitting 
model to analyze the experimental structure function. 


Second, “data smoothing” allows maximizing the software 
accuracy by making the new data shape as similar as possible 
to the learning data before choosing its best fitting equation. 
Third, “data clustering” predicts the more appropriate fitting 
equation of the experimental structure functions using a CNN. 
The number of possible fitting equations can be easily 
increased by adding more homogeneous-learning-classes in 
the learning phase. The approach is based on full-connected 
neural networks, which implies a voluminous calculation and 
very expensive computing time. Therefore, high-performance 
computing techniques using graphics processing unit (GPU) 
are highly recommended to obtain results in an optimal time. 
Once the software has selected the right model, an 
implemented Levenberg-Marquad non-linear least square 
fitting routine is used to extract the amplitudes and the decay 
times of the fluctuations. 
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Figure 2: Error between thermodiffusion experimental data 
and SF model as a function of the wave number using 
different fit equation, one exponential with constant and 
quadratic background, two exponentials with constant and 
quadratic background, represented by the gray, yellow, blue 
and green lines, respectively. The blue points represent the 
choice of the CNN program. 


Results 

In Fig. 2 the errors (averaged absolute value of the 
differences) between the experimental data and the fitting 
equations using four fit models (one exponential with 
constant or quadratic background and two exponentials with 
constant or quadratic background) are shown. The 
experimental data were obtained from Shadowgraph 
visualization of a thermodiffusion experiment on a 
polystyrene (PS)-toluene binary mixture. We used the 
concatenated angular average of bi-dimensional structure 
functions of three series of images acquired at 1Hz, 10Hz and 
100Hz. 


Our approach shows high performance on choosing the best- 
fit equations of the experimental structure function for each 
wave number. We notice that the software chooses the best 
fitting equation as compared to the minimum error for most 
of the wave numbers. 


Conclusions 

In this study, we present a new approach for automatically 
analyzing the experimental structure function using AI, 
particularly, artificial neural networks. Automatic selection of 
the fitting equation for each wave number greatly facilitates 
data analysis. It offers better performance in terms of result 
quality and analysis time, particularly in the presence of 
different physical phenomena and when analyzing a large 
amount of data. 
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Introduction 


The injection of CO2 into deep saline aquifers has been 
endorsed as one of the most promising pathways to reduce 
CO emissions in the atmosphere, as CO2 being considered as 
one of the major causes of the greenhouse effect (Bachu et al. 
2003). However, various studies on Carbon Capture, 
transport and geological Storage (CCS) have shown that 
injection of supercritical CO in a deep saline aquifer perturbs 
the physio-chemical equilibrium of the medium, leading to 
highly coupled Thermal-Hydrodynamical-Mechanical- 
Chemical processes (Kolditz et al. 2012). All these coupling 
effects are not yet well understood and formally implemented 
into the numerical models. In this framework, we are 
investigating the coupling between the thermal gradients and 
the reactive transport taking place in a diffusion cell, in order 
to model what happens in a saline aquifer. For this, an already 
developed diffusion cell is used (Ndjaka et al. 2022), where 
we put in contact two saline solutions and investigate the 
fluctuations happening during the mixing process. In this 
diffusion cell, the observations are performed in the direction 
parallel to the concentration gradients by means of 
shadowgraphy. Here, we present the results of the validation 
of the device, as well as the analysis method, on the free- 
diffusion of sodium chloride NaCl, sodium sulphate Na2SO. 
and calcium chloride CaCl: into water at 25 °C. By analysis 
of non-equilibrium fluctuations (NEFs), mass diffusion 
coefficients can be measured. A free diffusion experiment was 
carried out under homogeneous and _- stabilizing 
inhomogeneous temperature conditions and the impact of the 
temperature gradient on the mass diffusion process was 
investigated. After comparison of our results with literature, 
we were able to validate the experimental setup and 
methodology with the binary salt solutions. 


Experimental procedure 


The shadowgraph apparatus involves two main parts: a free 
diffusion cell, allowing to superimpose two layers of two 
different solutions and a precise differential temperature 
control through the liquid column, and the optical setup aimed 
at acquiring shadowgraph images of the density fluctuations 
(Croccolo et al. 2012). A light beam at wavelength (675+13) 
nm passes trought the sample in the vertical direction parallel 
to the imposed concentration gradients. The images acquired 
in real time by the camera (a CCD PIKE camera with a pixel 
side of l,j; = 7.32 um.) are saved on a workstation. In the 
absence of convective instabilities, we study the dynamics of 
fluctuations by calculating the Structure Function (SF) of the 
images by means of a differential dynamic Algorithm 
(Cerchiari et al. 2012). The expression of the normalized 
correlation function inthe SF fitting equation, that describes 
the time evolution of concentration fluctuations, has been 
modified in order to integrate the thickness of the sample. 


Series of N = 2500 images of 1024 x 1024 pix? have 
been recorded at frequency f = 7.8 Hz, every 5 minutes 
(after closing the valves) over the first 20 minutes of each 
experiment, and then every 10 minutes until the equilibrium 
state. 


Free-diffusion experiments with binary mixtures 


The cell filling protocol entails filling the diffusion cell 
completely with a less dense solution after which the bottom 
half of the cell is filled with a denser solution. This is achieved 
with the help of a syringe pump by filling the cell 
simultaneously with the two fluids from the bottom and the 
top inlets, while the remixed fluid is let out through the two 
outlets positioned at the mid-height of the cell. Once this is 
achieved, a relatively flat horizontal interface between the 
two solutions is formed. Then, inlet and outlet valves are 
simultaneously closed and the free-diffusion process starts. In 
order to validate setup, procedure and the adjustment model, 
we carried out a first series of measurements by 
superimposing two solutions of different concentrations but 
of the same salt. c is the mean molar concentration of 
dissolved salt between the two layers, which corresponds to 
the equilibrium concentration. Ac is the molar concentration 
difference between the bottom and the top layers. The average 
temperature of the sample was set at T = 25°C. Experiments 
were performed at atmospheric pressure. In Fig. 1, we plot the 
mean mass diffusion coefficient D as a function of the 
normalised time (i.e. time divided by the characteristic 
diffuvise time T,) for the three binary mixtures studied in this 
work, obtained after analysis of the NEFs decay times. 
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Figure 1: Mass diffusion coefficient D as a function of the 
normalized time for the free-diffusion experiments carried out at 
T = 25°C, blue and red circles for NaCl (c = 2.7 mol/L, Ac = 
2.0 mol/L ), blue squares for CaCl (c =0.46 mol/L, Ac 
0.8 mol/L) and blue triangles for Na2SOs (c = 0.5 mol/L, Ac 
0.7 mol/L). 


Literature values (Rard et al. 1979a, Rard et al. 1979b) are 
represented by the horizontal dashed lines in Fig. 1 for each 
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binary mixture of interest. The agreement is pretty good in the 
three cases. Red points in Fig. | correspond to a series of 
measurement carried out at the average temperature of T = 
25°C and a temperature difference between the top and the 
bottom of the column of AT = 20°C. The impact of the 
thermal gradient on the concentration NEFs was considered 
negligible. 


Conclusion 


We have studied by shadowgraphy the free-diffusion of the 
main salts representative of saline aquifers. An experimental 
setup has been developed and a first experimental procedure 
tested with saline solutions. We incorporated fluid column 
thickness effects into the fitting equation for the analysis of 
the measured structure functions. The procedure and the 
equation proved to be sufficient for the study of the free- 
diffusion process. 
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Introduction 


We propose a simple experiment that can be performed with 
household equipment, aimed at introducing students of all 
grades to the rich phenomenology of liquids stratified by 
diffusion, including processes such as levitation, Brunt— 
Vaisala oscillations and internal gravity waves. 


Stratification in fluids is often observed in the presence of 
density gradients controlled by diffusion and thermal 
diffusion and has a great relevance in many natural systems 
such as the atmosphere, the oceans, and the stars. 

A volume of fluid, or a body, with density in the same range 
as the stratification, can be in equilibrium at the appropriate 
height where the density coincides exactly with its own. 

A slight displacement from the equilibrium position causes 
the body to oscillate due to an elastic recovery force 
originating from the density umbalance. This oscillation, 
whose frequency depends on the density gradient and is 
known as Brunt—Vaisalé, can transmit its motion to the 
surrounding fluid and produce internal gravity waves. These 
phenomena are difficult to visualize as they only produce 
refractive index variations within the fluid. 

The physics involved can be quite easily understood and has 
many relevant links with diffusion, hydrostatic, physics of 
stratified fluids and internal gravity waves. For that reason, an 
experiment able to evidence these aspects has demonstrated 
to be a very good candidate for an engaging activity for 
students (Carpineti et al, 2021). 


In this new work, we propose a further experiment able to 
catch levitation, oscillations and gravity waves, but also the 
evolution of a density gradient through diffusion. The new 
experiment is a variant of the well-known Cartesian Diver 
experiment. In the traditional experiment, a hollow object is 
immersed in a sealed bottle containing a uniform liquid: when 
the applied pressure is varied the diver's density increases and 
the diver sinks, but when the original external pressure is 
restored, it goes back to the top of the liquid. 

In our experiment, the diver floats in a density stratified liquid 
obtained by pouring coarse kitchen salt inside a jar of water. 
The salt grains initially sink to the bottom of the jar, where 
they start to dissolve into water. The diffusion of the salt 
determines the formation of a gravitationally stable density 
gradient, which slowly evolves in time. By varying the 
applied pressure, the diver's density changes and it moves to 
a different height accordingly. In contrast to the original 
version of the experiment where the diver can only sink or 
float, in this case the diver can stop in a stable equilibrium 
position within the fluid, at the height where its density is 
matched. When a sudden pressure pulse is applied, the diver, 
pushed off its temporary equilibrium position, starts 
oscillating due to a restoring force. Moreover, by changing 
the pressure on the container, students can span different 


heights and, performing measurements at different times, they 
can map different density gradients and observe the evolution 
in time of the diffusion process with a single non-invasive 
experiment. 

It is also possible to make visible the internal gravity waves 
produced by the oscillating diver, by putting in suspension 
small fragments of a material of appropriate density, which 
localize in a fluid layer and oscillate when the gravity waves 
pass. With a simple experiment that students can project and 
realize by themselves with easy-to-find objects and by 
following all the steps, it is possible to let them experience 
how a simple diffusion process leads to complex phenomena 
of general interest. 

Thanks to the use of simple free educational programs, which 
allow recording the diver oscillations, plotting and fitting the 
data, students can perform quantitative analysis of the results, 
and therefore enhance their understanding of the physics 
issues. 


Figure 1 - A small glass bottle partially filled floats in 
equilibrium inside a density stratified liquid. Small plastic 
fragment with a known density are localized at a height where 
their density is matched and allow visualizing the internal 
gravity waves originated by the diver’s oscillations. 
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Introduction 


Chiral structures are frequently found in nature, on scales as 
diverse as those of elementary particles and living beings. 
These structures present mirror images known as 
enantiomers. In some situations, the ratios of these 
enantiomers are similar, but in others there is an excess of one 
over the other, resulting in chiral symmetry breaking. 
Explaining why this is so is a problem of great current interest 
(Hegstrom et al. 1990, Viedma 2005, Kondepudi et al. 2001). 


At equilibrium, the free energy AG,,, required for the 


formation of the possible enantiomers is practically the same 
AGrev 


so the corresponding probabilities p,~e *87 with kg the 
Boltzmann constant and T the temperature, are very similar, 
resulting in weak or none chiral symmetry breaking. 
Experiments have shown that the actuation of external factors 
such as temperature gradients (Ribo et al. 2013) and in 
general external forces driving the system out of equilibrium 
can lead to a significant disproportion in the concentrations of 
the homochiral structures (enantiomers) or even prevent their 
formation (Buhse et al. 2021, Szurgot 2012). 

In this paper, we demonstrate by modelling and experiments 
that the energy required to generate an enantiomeric crystal 
excess must be equal to the energy dissipated in the 
irreversible processes involved in crystallisation: heat 
exchange and salt phase change. This energy causes an 
increase in the free energy barrier of one of the enantiomeric 
crystals and, consequently, a decrease in the enantiomer 
formation rate which explains the enantiomeric excess. In the 
experiments performed, we have measured for the first time 
the enantiomeric excess of NaClO3 along the crystallisation 
process of this substance, for different non-equilibrium 
conditions. We have found that the enantiomeric excess 
becomes more important when the actuating forces, 
temperature and activity differences, increase and therefore 
when energy dissipation, measured in terms of the entropy 
production rate, is larger. 


Chiral symmetry breaking 
production 

The existence of an enantiomeric excess has been attributed 
to the disparity in the activation energies of both enantiomers 
which causes different formation rates (Rib6 et al. 2017, 
Kondepudi et al. 2001). In Fig. 1, we illustrate the different 
stages of crystallisation in which the two driving forces, 
temperature difference AT between the system and its 
surrounding and activity difference Az between the solid 
and liquid phases reach different values. 


induced by entropy 


Az=0 Az>0 Az=0 

AT = 45°C AT = 30°C AT =0°T 
Figure 1: Enantiomer formation process. The system is initially an 
under-saturated mixture of salt, denoted by small black semicircles, 
and glycerol. Temperature decrease leading to stage a) characterised 
by stronger interaction among salt compounds (yellow and blue 
lines) thus reaching saturation. Small nuclei of both enantiomers are 
formed (L- and D-nuclei represented by gold and blue agglomerates 
respectively) in the process of phase change in which one of the 
enantiomers is dominant, as represented in stage b). Finally, crystals 
emerge and the system reaches a state of chemical and thermal 
equilibrium c) with an enantiomeric excess. 


Theory 
Heat transfer 
The heat that the system exchanges with the environment 
which is at constant temperature T, causes the temperature 
varies according to the equation 

OT 0°T (1) 


U 
i agg Ce Baye 


Where c, is the specific heat of the mixture at constant 
pressure, U the convective heat transfer coefficient, cg the 
total molar concentration, d the test tube diameter, and k 
the thermal conductivity. 


Crystallisation kinetics 
The evolution equation of the average L-crystal size a, is 


(2) 


2 
da, s Dlg (52) : pa, (=) (Ar)? 
dt Kgl \Oa, Kp \ OT 

With ¢, the potential energy to form a crystal of size a,, 
D, the diffusivity in the | space (Reguera et al. 2005), lo 
the minimum size of the nucleus, n, the amount of salt in L- 
configuration, k,, the thermal expansion coefficient for the 
solid salt and % the coupling constant. An analogous 
evolution equation can be written for D-crystals. 


Results 

We obtain Fig. 2 for one initial and cooling condition. We 
observe that heat transfer and nuclei formation are the most 
dissipative processes and that the supersaturation gradient is 
the main driving force for crystal formation (Szurgot 2012b), 
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as it causes the most dissipation. This conclusion is reached 
in all the cases studied. 
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Figure 2: Energy dissipation rate corresponding to the different 
irreversible processes occurring in crystal formation as a function of 
time. 
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Figure 3: Energy dissipation rate corresponding to the different 
irreversible processes occurring in crystal formation as a function of 
time 

In Fig. 3(a-b), we represent the energy dissipation rate of: i) 
heat transfer with the environment due to a temperature 
difference; ii) solid nuclei formation driven by the over- 
saturation force; iii) crystal growth caused by surface and 
volume energy differences. Fig. 3(a) shows the presence of 
two peaks and one minimum whereas in Fig. 3(b) there is only 
one peak for each salt/solvent concentration ratio. We also 
observe that for enhanced cooling the peaks are slightly 
higher than those for natural cooling which is a consequence 
of the fact that in the former case, the energy barrier is higher 
and consequently the barrier crossing rate diminishes which 
leads to an increase of the enantiomer population difference, 
as shown in Fig. 3(c-d). The final enantiomeric excess 
percentage is higher in the case of enhanced cooling, as the 
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energy dissipation rate is initially higher. Therefore, the early 
excess of one of the enantiomers, with a lower redissolution 
rate, causes the enantiomeric crystals to persist over time (A. 
Arango-Restrepo et al. 2023). 


Conclusions 

The minimum values of the energy dissipation are found 
when the dissipation due to phase change decreases and heat 
exchange becomes the most dissipative process. This is 
connected with local maximum of the enantiomeric excess 
which means that temperature gradients might amplify the 
symmetry breaking. 

The enantiomeric excess depends on the cooling rate and 
temperature gradients. The path in which the cooling process 
takes place, inducing supersaturation of the mixture, affects 
the intermediate states in the nucleation, thus leading to 
chiral symmetry breaking. 


Acknowledgements 

A.A-R. and J.M.R. are grateful for financial support of 
MICIU (Spanish Government) under grant No. PID2021- 
126570NB-I00. O.A. recognises the financial support of 
MICIU (Spanish Government under grants No. RYC2018- 
024997-I and RTI2018-098410-J-I00. 


References 

R. A. Hegstrom and D. K. Kondepudi, The handedness of the 
universe, Scientific American, 262, 108-115, (1990). 

C. Viedma, Chiral symmetry breaking during crystallization: 
complete chiral purity induced by nonlinear autocatalysis and 
recycling, Phys. Rev. Lett., 94, 065504, (2005). 

D. K. Kondepudi and K. Asakura, Chiral autocatalysis, 
spontaneous symmetry breaking, and stochastic behavior, 
Accounts of Chemical Research, 34, 946-954, (2001). 

J. M. Ribd, J. Crusats, Z. El-Hachemi, A. Moyano, C. Blanco 
and D. Hochberg, Spontaneous Mirror Symmetry Breaking in 
the Limited Enantioselective Autocatalysis Model: Abyssal 
Hydrothermal Vents as Scenario for the Emergence of 
Chirality in Prebiotic Chemistry, Astrobiology, 13, 132-142, 
(2013). 

T. Buhse, J.-M. Cruz, M. E. Noble-Teran, D. Hochberg, J. M. 


Ribé, J. Crusats and J.-C. Micheau, Spontaneous 
deracemizations, Chemical Reviews, 121, 2147-2229, 
(2021). 


M. Szurgot, in Crystallization and Materials Science of 
Modern Artificial and Natural Crystals, InTech, Rijeka, 
Croatia, (2012). Parity Violation in Unstirred Crystallization 
from Achiral Solutions. 

J. M. Rib6, D. Hochberg, J. Crusats, Z. El-Hachemi and A. 
Moyano, Spontaneous mirror symmetry breaking and origin 
of biological homochirality, Journal of The Royal Society 
Interface, 14, 20170699, (2017). 

D. Reguera, J. M. Rubi and J. M. G. Vilar, Mesoscopic non- 
equilibrium thermodynamics, The Journal of Physical 
Chemistry B, 109, 21502—21515, (2005). 

M. Szurgot, Chiral symmetry breaking in unstirred 
crystallization, Crystal Research and Technology, 47, 109-— 
114, (2012). 

A. Arango-Restrepo, O. Arteaga, D. Barragan, J.M. Rubi, 
Chiral symmetry breaking induced by energy dissipation, 
PCCP, submitted (2023). 


79 


15" International Meeting on Thermodiffusion (IMT15) 


Tarragona (Spain), May 2023 


Simultaneous appearance of fingers and overstable instability in isothermal ternary systems 


Berin Seta', Jon Spangenberg!, Fina Gavalda2, M. Mounir Bou-Ali?, Xavier Ruiz’, Valentina Shevtsova?“ 


‘Department of Civil and Mechanical Engineering, Technical University of Denmark, Lyngby, Denmark 
"Department of Physical and Inorganic Chemistry, Universitat Rovira I Virgili, Tarragona, Spain 
Fluid Mechanics Group, Mondragon University, Arrasate/Mondragon, Spain 
‘IKERBASQUE, Basque Foundation of Science, Bilbao, Spain 


Introduction 


Buoyancy can be generated by multiple driving forces 
simultaneously, as in the cases of temperature and 
concentration gradients or gradients of two or more species in 
multicomponent mixtures. The convective movements due to 
density inversion caused by different rates of diffusion are 
called double-diffusive convection. Double-diffusive 
convection appears in many natural processes, most notably 
in form of salt fingers causing water mixing in oceans (Stem 
1960). However, the convective motion caused by different 
diffusion rates of species disturbs experimental 
measurements in many different setups, most notably the ones 
with a sharp interface between solutions, like Sliding 
Symmetric Tubes (SST), Counter-Flow Cell (CFC) or Open 
Capillary Tubes (OECT) (Seta et al. 2019). 

In ternary mixtures, it is well established that two different 
types of instability can occur: fingers-like, where “fingers” of 
one mixture are transported into another mixture and 
overstable type of instability, where the sharp interface is 
maintained through experiments, but convective structures 
are formed on the edges of interface. Static limits forthe onset 
of previously mentioned instabilities can be obtained either 
by taking derivatives of the analytical solution of the diffusion 
transport equation in the ternary mixture or applying linear 
stability analysis (Seta et al. 2019, Vitagliano et al. 1992). 
Nevertheless, it must be mentioned that the inversion density 
pockets are sufficient, but not required for the appearance of 
finger-like instability. 

The experiments conducted in the aforementioned setups are 
designed to maintain a small concentration difference 
between the top and bottom mixtures. This is done to ensure 
that the constant diffusion coefficient assumption is valid. 
However, some experiments have been carried out with a 
concentration difference of up to 0.1 (Mialdun et al. 2013, 
Larranaga et al. 2014). This difference is significant enough 
to result in distinct diffusion coefficients across the 
concentration gradient, as demonstrated by the ternary 
diagram in Figure 1 forthe Toluene (Tol) — Methanol (Meth) 
— Cyclohexane (Cyl) mixture in the case of Diz. 

By using a set of constant diffusion coefficients, stability 
graphs can be created based on the initial concentration 
difference between independent components toluene and 
methanol. These graphs can indicate the initial conditions that 
will result in convective motion. However, if we account for 
the concentration-dependent nature of diffusion coefficients, 
it is possible that the top and bottom mixtures will have 
different coefficients. In this case, the stability graph for each 
mixture may differ, leading to uncertainty in predicting the 
occurrence of instability during the experiment. 
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Figure 1: Concentration-dependent variation of Diz. Red 
rectangular shows a region of interest, where, for a relatively small 
variation of concentration, cross-diffusion term Dj2 shows big 
differences. 


Methodology 


The system is initialized as two superimposed layers of 
ternary mixture with slightly different concentration and 
hence different diffusion coefficients. The Navier-Stokes 
equations, along with the transport equations for independent 
components wi (toluene) and w2 (methanol), are solved using 
the open-source software OpenFOAM. The mixture is 
assumed to be incompressible, with the exception of the 
Boussinesq approximation, where the density depends on the 
mass compositions. The governing equations can be 
expressed as follows: 


tu Vu=——+Vpt+vAut gp’ (1) 
ot Po 

Mt +u-Vw, =V-(yVW,)+V-DpzVw,) 2) 
M2 + u- Vy = V- (DV) +V-(DzVw2) 3) 


Where wu is velocity vector, Di; are concentration-dependent 
diffusion coefficients, p) is mean density, p* = p)(1 — 
By1 Wi — Wip) — By2 (wz —W29)) is mixture density 
dependent on concentration of species, B,,, and B,,. are mass 
expansion coefficients of toluene and methanoland g is 
gravity vector. The stability graphs are obtained from limits 
given in (Seta et al. 2019) for constant diffusion coefficients. 
Simulations were conducted in Hele-Shaw cell. 


Results 


The initial difference in diffusion coefficients between the top 
and bottom solutions in the Hele-Shaw cell results in an 
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asymmetric solution (Bratsun et al. 2022). However, the 
influence of cross-diffusion in this work was neglected, which 
is crucial mechanism to obtain simultaneously two different 
types of instability. Furthermore, previous research (Bratsun 
et al. 2015) suggested that another coupling mechanism, such 
as a reaction, was necessary to produce two different 
instabilities simultaneously in different parts of the diffusion 
front. In this study, we demonstrate that if the diffusion 
coefficients are such that the stability pattern undergoes a 
complete change in the top and bottom solutions, meaning 
that the regions of overstable and fingers switch places, it is 
possible to observe two _ instabilities occurring 
simultaneously. For example, if we initialize system with 
W1,bot=0.67 and W2,bor=0.015 and Wi top=0.65 and W2.top=0.065, 
maps of stability will look like in Fig. 2. Based on the work 
of Grossmann (Grossmann et al. 2009), polynomial equations 
are used to fit the diffusion coefficients as a function of 
concentration. The regions of overstable and fingers have 
evidently shifted, resulting in distinct movements in different 
areas of the domain. 


a) Top solution b) Bottom solution 
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‘1 -0.05 0 0.05 01 0.1 0.05 0005 
W2,t-W2.b W2,t-W2,b 

Figure 2: a) Stability map for top solution, b) Stability map for 

bottom solution. 


Once the conditions for different types of instability are 
reached at the top (fingers) and bottom (overstability), the 
resulting fluid motion differs, as shown in Figure 3. The 
resulting structure is asymmetric, with small fingers forming 
in the upper part and large structures forming in the bottom, 
far from the center and hindering movement, which helps to 
maintain a sharp interface. 
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829.8 
t= 600s 
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overstable 
. 830.1 
pattern in the 
bottom with a 830 
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moved to the 829.8 
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Figure 3: Formation of two different instabilities in a ternary 
mixture with diffusion coefficients depending on the composition. 


Formation of 
fingers in the 
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As seen in the later stages of the process (t=1250s), the entire 
interface moves downward, as the upper part is mixed due to 
convective motion of the fingers, while the bottom remains 
blocked. Furthermore, the onset of instabilities is time- 
dependent and differs at early stages, suggesting that one 
instability may evolve into another once the field is 
sufficiently mixed. 


Conclusions 


This study presents, for the first time in literature, the 
simultaneous appearance of two types of instability in an 
isothermal ternary system. The primary mechanism behind 
the asymmetric solution is the concentration-dependence of 
the diffusion coefficients. This concentration-dependence 
allows for the emergence of various patterns, such as finger- 
overstable, finger-stable, stable-overstable, and others. 
Additionally, the onset of these patterns may occur at 
different times and be transient in nature, causing certain 
instabilities to emerge and disappear during the course of the 
experiment. 
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Introduction 


The aim of the present paper is to investigate the linear 
stability of the Poiseuille flo w of binary fluids by taking into 


account both the effects of viscous heating and 
thermodiffusion. Both boundaries are considered 
impermeab le. The upper boundary is assumed 
isothermal while the lower one is adiabatic. 


Therefore, this work may be seen as an extension to binary 
fluids of the linear stability analysis conducted by [Barletta 
et al. 2011] for a one-component fluid. As no external 
temperature or concentration difference is 
imposed on the layer, the cause of thermal 
instability is the flow rate through the 
volumetric heating induced by the’ viscous 
dissipation and the Soret effect inherent to 
binary mixtures. 


Physical problem and the basic state 

We consider a binary fluid mixture contained 
between two horizontal plates of _ infinite 
extension in the horizontal directions xand y, 
separated by a height H in the vertical 
direction z, so that the lower wall is in z=0 
and the top wall in z=A. A constant pressure 
gradient is imposed in the x -direction, 
generating a Poiseuille velocity profile. The 


thermal boundary conditions are such that the 
lower wall is adiabatic while the upper wall is 
isothermal. No slip for the velocity field and 


no mass flux conditions are imposed on _ the 


horizontal walls (Fig. 1). 


(j.7= 0) 


(Homogeneous fluid) 
mixture Pix 
Umax —pP 


(Pin > Pout) 


Figure 1: Geometry problem 


The basic state solutions may be written in dimensionless 


form as: 


u, = Re Pr f(z )es3 T, =A g(z), C,=-y i, 
The functions f(z) and g(z) are defined as follows: 


f(z)=6z(1-2z), g(2) = (2-32 +424 —2z") 


where Re is the Reynolds number, Pr is the Prandtl 
number, A is the viscous dissipation Rayleigh number 
defined as A=Ge Re’ Pr* with Ge the Gebhart 
number, a parameter that measures the intensity of the 
viscous dissipation and the parameter yw is the 
separation ratio which represents the ratio 
between the mass contribution and the 
temperature contribution to buoyancy forces. 

The above basic solutions state that in the presence of an 
imposed Poiseuille flow, a negative basic temperature 
nonlinear vertical gradient exists due to the viscous 
dissipation contribution. A basic vertical nonlinear 
concentration stratification is then generated by the Soret 
effect. 


Linear stability of the basic state 

A linear stability analysis of the above basic 
solutions has been carried out in [Ali Amar et 
al. 2022] by solving the eigenvalue problem 
numerically using the shooting method for 
arbitrarily oriented convective rolls. The 
results indicated that the most unstable 
disturbances are longitudinal rolls 
independently of the Soret effect. In addition 
to Lewis number (Ze) and Prandtl! number Pr , 
the onset of convection in the’ form of 
longitudinal rolls depends on the viscous 
dissipationRayleigh number A. 


Results for yw > 0 


It is interesting to note that both the basic temperature and 
the basic concentration gradients are negative, meaning that 
the system may be unstable under some conditions. We 
present in Fig. 2(a) the linear stability diagram in the 


(w, A) plane. 
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Figure 2: Stability diagram for yw > 0 


The stability diagram depends on the values of the Lewis 
and Prandtl numbers. The plot in Fig. 2(a) was obtained for 


Pr =10 and different values of Lewis number Le = 50, 
Le =100 and Le = 200, values which are typical for 
liquid mixtures. The same figure shows the stability diagram 
for Pr=1and two values of Lewis number Le =1 and 
Le = 2, representing gas mixtures. As can be seen in this 
figure, the instability appears ata lower A. than in the case 
of a mono-constituent fluid (yw =0) for a liquid mixture 


and fora gas mixture. The difference is that the destabilizing 
is more important fora liquid mixture than fora gas mixture. 
Moreover, this figure shows a drastic drop of the critical 


threshold A, as the Lewis number increases. Numerical 
results indicate that instability appears with zero frequency, 
indicating that the emerging longitudinal rolls are stationary. 
As shown in Fig. 2(b), for a fixed y #0 the instability 
appears with a wave number k, that deviates weakly from 
its value for y = (in the case of a gas mixture, while it 
decreases strongly for a liquid mixture as yw increases. We 


remark that the decrease in the wave number is drastic for a 


liquid mixture with a high Lewis number. In that case, there 
exists a particular value of the separation ratio Wig», Such 
that a finite critical wave number is observed for 


W < Wirono and zero wavenumber for W > Winono - 


Results for Wr 0 


For W ~ Q, the basic temperature and concentration fields 
are of the same sign and the densest constituent migrates 


toward the hot wall, that is, the lower wall, stabilizing the 
system. Figure 3(a) illustrates the evolution of the threshold 


A, at the onset of convection as a function of y for 
different values of Le and Pr. For fixed Le and Pr, 


two types of instability may exist depending on the value of 
y . On one hand, the instability appears with a non-zero 


wave number and non-zero frequency if yw is below a 


particular value Wc; . Thus, the instability mechanism is an 


oscillatory periodic Hopf bifurcation, leading to symmetry 
degenerating left and right traveling wave patternsin the y 


direction. At the onset of convection, their corresponding 
threshold A. and frequency @, are plotted as functions of 
separation ratio by the dashed lines of Fig. 3. The dashed 
lines of Fig. 3(a) show that A, increases very weakly with 
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Figure 2: Critical values at the onset of stationary longitudinal rolls 
(solid line) and oscillatory longitudinal rolls (dashed line): (a) A,, 
(b) @ vs yw for Le=100, and Pr=10 (square symbols), 


Le=2 and Pr=1 (dot symbols) and Le = Pr =1 (diamond 
symbols). 


We now restate the results obtained in terms of the critical 
Reynolds number that can be easily evaluated from the 


relationship A = Ge Re’ Pr’, namely, 


Reowe ae 
Pr \ Ge 


We note that Re, is proportional to the inverse of the 


Prandtl number, which means that the more we increase the 
Prandtl number, the more the critical Reynolds number 
decreases. This may allow us to observe the effect of viscous 
dissipation on experimentally affordable Reynolds number 
for highly viscous fluids. Figure 4 presents the stability 
diagram in the (yw, Re) plane for fixed Pr and Ge. We 
observe in Fig. 4 that the larger the positive separation ratio, 
the smaller the critical Reynolds number needed to trigger 
stationary convection, which confirms once again the 
destabilizing role of the Soret effect with positive separation 
ratios. Inversely for negative separation ratios, the Soret 
effect hasa stabilizing influence. 


w=0, Le=2 


BET 


' 
6 

\ 
2 
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Figure 4: Critical value Re. vs yw for Le=2, 100 , 


Pr=10 and Ge=10~ at the onset of longitudinal rolls: 
stationary mode (solid line) and oscillatory mode (dashed line) 


Conclusions 

The linear stability analysis is performed for the convective 
instability induced by viscous dissipation for a Poiseuille 
flow of a binary fluid mixture by taking into account the 
Soret effect. It is found that 

(i) For positive y , the bifurcation is stationary, and the 
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Soret effect precipitates the appearance of stationary 
longitudinal rolls. Moreover, a drastic drop of the critical 
wave number is observed for liquids such that it becomes 
zero if y exceeds a particular value W,,5,.- This means that 
longitudinal rolls have an infinite wavelength or are 
structured as a mono-cellular flow in a real channel bounded 
in the y direction; 

(i) For negative y , the system exhibits a Hopf bifurcation, 


if wy < Wcr . Otherwise, the bifurcation is stationary. 
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Introduction 


In recent years, the interest and study in the phenomenon of 
thermodiffusion that describes the separation of species under 
a temperature gradient has grown significantly due to its 
applicability in different sectors (W. KGhler et al. 2016). The 
representative magnitude of this phenomenon is the Soret 
coefficient (St) which in a binary mixture is the fraction 
between the thermodiffusion (Dr) and diffusion (D) 
coefficient. When Sr is positive, the denser component of the 
binary mixture moves to the cold region. Otherwise, in a less 
common scenario, the denser component goes to the hot part. 
Regarding the experimental process, the thermogravitational 
column (TGC) technique was one of the first methods 
applying the Dr effect under thermogravitational conditions 
to separate uranium isotopes (K. Clusius et al. 1939). The 
analytical solution to the problem was presented by W. H. 
Furry, R. Clark Jones and L. Onsager and is known as Fury- 
Jones-Onsanger (FJO) theory (W. H. Furry et al. 1939), who 
succeeded in associating the thermodiffusion coefficient Dr 
of a binary mixture from stationary separation measurements 
Ac, between two points AL, using the following Eqn. (1). 


Vv 
Ac = 504 eo — Co) (1) 


where v is the kinematic viscosity of the sample, o the 
thermal expansion coefficient, g is the gravity force, cg the 
initial mass fraction and L, the gap of the column. 


Binary mixtures of positive Soret coefficient have been 
studied in detail both numerically and experimentally using 
the cartesian configuration TGC technique (Seta, B. et al. 
2020). While in the case of the negative ones, the context is 
more complex. There is an experimental work where an 
adverse density gradient was stabilised in a TGC of 
cylindrical configuration (M. M. Bou-Ali et al. 1999), 
including the determination of negative Dr coefficients in 
different binary mixtures (M. M. Bou-Ali et al. 2000). Later, 
(Kozlova, S. et al. 2020) analysed the thermohydrodynamic 
stability problem under the same conditions as (M. M. Bou- 
Ali et al. 1999) stabilising an adverse density gradient by 
considering a larger gap in the cyclindrical configuration. 
This problem was also analysed for multicomponent mixtures 
(Ryzhkoy, I. I. et al. 2009), pointing out that the long-wave 
perturbations were less invasive than for the binary case. 
Moreover, (Zebib, A et al. 2008) concluded that transversal 
waves (in third dimension) were more dangerous than 
longitudinal ones. Recently (Berin Seta et al. 2020) 


investigated the stability analysis under thermogravitational 
effect concluding that the third dimension (not considered in 
the FJO theory) has a significant effect on the stability of an 
adverse density gradient in the cartesian configuration TGC. 
In this work, a numerical analysis was carried out to 
determine the best geometrical conditions of the cartesian 
TGC configuration to stabilise a negative binary mixture. A 
new TGC installation was designed and fabricated to verify 
and validate the numerical results. 


Numerical model 


As a Starting point, the thermogravitational microcolumn 
(micro-TGC) existing in our laboratory was used. This has a 
gap of 0.509 mm (Lx), a width of 3 mm (Ly) and a height of 
30 mm (L,). Ansys Fluent 2022 R1 software was used 
considering the species transport model to simulate the 
species separation of a binary sample in the micro-TGC. The 
boundary conditions were defined in such way that a 
horizontal temperature gradient was applied between two 
parallel walls, while rest of the walls were considered 
adiabatic. As for the properties of the mixture, it was defined 
that the density of the sample changes due to temperature and 
concentration considering the Boussinesq approximation Eqn. 
(2) where Pg is the initial density of the sample, Tp is the 
mean temperature and B is the mass expansion coefficient. 


P = Po[1 — a(T — To) + B(c — co)] (2) 


For the validation of the numerical model, the mixture 
tetrahydronaphthalene-dodecene (THN-Ci2) at 50 % in the 
mass fraction and at 25 °C was considered as a reference 
(Platten, J. K. et al. 2003). The variation of the concentration 
of the denser component (THN) along the micro-TGC was 
analysed for AL, = 25 mm. As can be seen (Table 1), the 
deviation between literature, experimental and numerical 
results do not exceed 5 %. 


Table 1 : Stationary Ac values sample THN-Ci2 50 % at 25°C. 


Acret= 0.0576 AcExp = 0.0583 Acnum = 0.0579 


Once numerical model was validated, the thermogravitational 
behaviour was analysed in the case of a mixture with a 
negative Soret coefficient, toluene-cyclohexane (tol-Ch) 
50 % mass fraction at 25°C. Different widths of the column 
were analysed, Ly = 3, 1.75, 1.5 and 1 mm, always 
maintaining the height and the gap of the micro-TGC as 
constant (L, = 0.509 mm and L, = 30 mm). The numerical 
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results provided that the effect of the third dimension 
stabilizes the system for the widths of 1mm and 1.5 mm. (see 
Figure 1). The simulation used the value of the Soret effect 
given in the literature (Lapeira et al., 2017). It was noted that 
separation in a column with a width of 1 mm had a higher 
value when compared to FJO theory (literature), so the 
cartesian configuration of a 1.5 mm window was considered 
the optimum for that mixture. 


x10% 


'— Literature 


20 40 60 80 100 120 
t/ min 


Figure 1: Separation of the sample Tol-Ch 50 % at 25°C. Numerical 
results for different geometrical configurations. 


Experimental setup and preliminary results 


Based on the numerical results, a new micro-TGC was 
designed, fabricated, and adapted to the optical installation in 
the laboratory, see Figure 2 a). The blue laser (473 nm) was 
used in the experiments, which has shorter wavelength when 
compared to the red one (633 nm) and allows us to obtain a 
higher phase separation. An optical sensor based on a Mach- 
Zehnder interferometer was used to measure the change in 
concentration between two selected points. The experimental 
results shown in figure 2 b), confirm that, for the first time, an 
adverse density gradient was experimentally stabilised in a 
binary mixture validating the results shown in the numerical 
section corresponding to the mixture tol-Ch 50 % at 25°C. 


” 
y 
: 


: 
pesooboad scene 


200 400 600 800 1000 
t/ min 


Figure 2: New fabricated micro-TGC a) and experimental results, 
separation of the mixture tol-Ch 50 % (stable) at 25°C b). 


Conclusions 


In this work, the effect of the third dimension on the stability 
of a separation with an adverse density gradient under 
thermogravitational effect in a cartesian configuration was 
analysed. For the first time, it has been numerically and 
experimentally demonstrated that an optimal choice of TGC 
width allow to measure a binary system with a negative Soret 
effect otherwise unstable. To do so, a new micro-TGC was 
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designed, fabricated and set up by changing the column width, 
also known as third dimension. 
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Introduction 


Carbon dioxide (CO2) long-term storage has been identified 
as one of the most relevant strategies to achieve the 2015 Paris 
Agreement objectives (Rogelj et al., 2021). Although the 
main physico-chemical processes have been well 
characterised (Emami-Meybodi et al., 2015; De et al., 2022), 
several open questions remain regarding the 3D convective 
dissolution of CO in brine. This is mainly because of the lack 
of experimental data close to process-relevant conditions. 

As a first step towards the understanding of convective 
dissolution of COz in brine, we investigated the convective 
patterns in a 3D system in free medium (Fruton et al., 2023). 
More recently, we started preliminary study to observe the 
same phenomena in a model transparent porous medium. In 
this context, we developed an experimental setup using the 
shadowgraph method for the dynamic spreading of 
convective plumes in saturated transparent porous media, 
whose permeability ranged between 6 and 60 x 10°!° m?. The 
first studies were performed close to the refractive index 
matching (RIM) conditions without fluorescent dyes using 
fluid mixtures with density gradients similar to those involved 
during the dissolution of CO: in brine. 


Experimental Methodology 


Here, the optical technique is shadowgraphy. The optical 
apparatus consists of a super-luminous diode (SLD) and a 
high-resolution s-CMOS camera from Hamamatsu. The light 
beam is passing through the sample cell which is a 3500 pL 
precision cuvette (Thorlabs). It contains a 3D model porous 
medium filled with a matching solution of 1-hexanol and 
toluene as shown in figure 1. The porous medium consisted 
in borosilicate glass beads (Sigma-Aldrich) of 1 or 3 mm 
diameter. 

A denser fluid was injected through a syringe pump by 
allowing a drop of 16.7 wL volume to fall at the liquid-gas 
interface. The drop is prepared to have a density higher than 
the solution saturating the porous medium. Three density 
differences (Ap=13.6, 7.7, and 4.9 kg/m*) were tested to 
investigate the convective process through the transparent 
porous medium. 

During the experiment, we ensured that the drop was injected 
very slowly and that the injection needle was placed at the 
centre of the cuvette. Pure fluid density and viscosity data 
were obtained from the supplier and compared to literature 
data (Wu et al., 2019; Almasi, 2021). Series of images were 
acquired at a frequency of 20 Hz and the experiments 
typically lasted for about 5 minutes. We defined Z=0 for each 


experiment at the liquid-gas interface as shown in Fig. 1.To 
test the reliability of our setup, each experiment was repeated 
at least three times. 


Drop fluid 


Plumes 1 mm beads 
direction 


Outside FOV* 


3 mm beads 


*FOV - Field of view of camera 
Figure 1: Schematic representation of the cuvette filled with a 
porous medium saturated with the RIM solution. The free liquid-gas 
interface is chosen as the origin of the vertical axis. On the right- 
hand side, the shadowgraph images show the position of the porous 
medium and liquid-gas interface. 


Results and Discussions 


For each density difference, we analysed the convective front 
velocity (CFV) defined as the velocity of the first plume 
travelling downwards to the bottom of the sample cell. We 
computed the variance of the lines of each image difference 
at different times after injection. Using a threshold variance, 
we were able to determine the position of the first plume, Zo, 
see figure 2. 

For both porous media, the CFV increases as Ap increases. 
While average CFV doubles when Ap increases from 7.7 
kg/m? to 13.6 kg/m? for the 1 mm beads (0.29 mm/s to 0.59 
mm/sec), the speed increase was less than 20% for 3 mm 
beads (3.43 mm/s to 4.18 mm/s). For the experimental 
dissolution of COz in brine at a density difference of 0.29 
kg/m3, the maximum vertical velocity of the leading fingers 
in a Hele-Shaw cell is approximately 0.029 mm/s (Li, Zhong 
and Jiang, 2023). 

Our results show that when a large Ap is achieved (1.5%), the 
plumes travel faster, but when Ap is small (0.9% and 0.6%), 
the spreading velocity is affected by the tortuosity pathway of 
the porous media. The density difference attained during 
convection substantially influences the CFV. 
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Figure 2: CFV analysis for the 1 mm porous medium and a density difference of 4.9 kg/m°. Top left: image difference at t=10 s. Top right: 
variance of each line of the image difference. The intersection of the line variance curve (blue) and a threshold variance (black line) determines 
the position of the plumes, Zo, at a given time. Bottom: plume position (blue), Zo, and velocity profile (orange). The average CFV of this 


system's early plumes was 0.27 mm/s. 


Conclusions 


In this study, we developed an experimental methodology to 
study density-driven convection in a model transparent 
porous medium where we could experimentally determine the 
convective front velocity for three values of Ap = (4.9, 7.7 and 
13.6 kg/m*) and two permeabilities of 6 x 10° and 6.6 x 10°'° 
m’. The CFV analysis shows consistency and repeatability. 
We show that when a high Ap is achieved, the early plume 
velocity is significantly affected by the density difference, but 
at Ap lower than | %, the early spreading velocity of the 
plume considers other factors of the porous media (e.g. pore 
spaces). These results can be used as input parameters for 
simulations of convective mixing in porous media. The 
relationship between the laboratory experimental data and its 
implications for CO sequestration will be explored in the 
near future. 
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The Marangoniconvection in a heated layer of a surfactant 
solution in the presence of Soret effect has been a subject of 
research during past decades; forreview, see (Shklyaev et al. 
2017). Typically, it is assumed that the concentration of the 
surfactant is below the critical micelle concentration. It is 
known that the influence of the micelle development on the 
Marangoni stresses is significant (Edmonstone et al. 2006, 
Craster et al. 2009). 


We consider the development of the Marangoniconvection n 
a heated layerofa surfactant. The Soret effect, the kinetics of 
adsorption and desorption of the surfactant on the liquid 
surface, and the kinetics of micelle formation are taken into 
account. The polydispersity of micelles is neglected. The 
main attention is payed to the cases where a longwave 
Marangoniinstability is developed (moderate Galileo number 
or small Biot number). The longwave expansions are used for 
the derivation of strongly nonlinear equations for active 
variables. The linear stability theory is developed for 


monotonic and oscillatory instability modes, and the pattem 
selection is discussed. 
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Introduction 

By applying a temperature gradient to a binary fluid mixture, 
a concentration gradient is generated by the Soret effect. In a 
non-associative system, like the binary mixture of |-hexanol 
and carbon dioxide (COz), the Soret coefficient is negative 
which can lead to hydrodynamic instabilities in the form of 
convective rolls. Before convection sets-in, non-equilibrium 
fluctuations (NEFs) spontaneously emerge in the presence of 
a density gradient and can be additionally amplified by 
buoyancy, thus originating convection (Giavazzi and Vailati 
2009). The Shadowgraph method allows to visualize and 
analyze both NEFs and convective patterns (Croccolo and 
Brogioli 2011). The signal stemming from convective 
patterns, however, overwelhms the one from NEFs. Thus, 
studying this process under reduced gravity conditions can 
provide a deeper understanding of the behaviour of NEFs in 
conditions close to the onset of convection. 

In this work, we present the set-up specifically developed for 
the 61" CNES parabolic-flight campaign VP161. Then, we 
describe the analysis method and the results obtained during 
thermodiffusion experiments at different gravity levels (1g, 
1.8g and 0g). The samples are mixtures of 1-hexanol and COz 
at different molar concentrations. Positive and negative 
temperature gradients are applied to the system to study both 
stable and unstable conditions. The influence of the different 
gravity levels on convection is investigated. 


Materials and Method 

The sample under study is prepared in a vapour-liquid- 
equilibrium (VLE) vessel by mixing CO: with 1-hexanol at a 
given pressure. The vessel is strongly shaked to reach a VLE 
quickly with the concentration of CO: in the 1-hexanol phase 
determined by the thermodynamic conditions of temperature 
and pressure. 

Thereafter, the liquid phase of the mixture is transferred to a 
sample cell that mainly consists of a stainess-steel annulus 
sandwiched between two  thermally-controled sapphire 
windows. The thermal gradient is applied by two Peltier 
elements featuring a 13 mm-diameter hole for the optical 
access. The cell is placed in a Shadowgraph apparatus 
consisting in a beam, provided by a superluminescent diode 
(SLD) and collimated by a lense, passing through the sample 
cell towards a high-speed sCMOS camera. The latter can 
record intensity maps of 2048x2048 pixels at a maximum 
frequency of 100 Hz, thus providing information about the 
variations of the refractive index and, thus, of the density 
variations inside the cell. 

During an experiment, we first apply a thermal gradient, and, 
once the macroscopic thermal equilibrium is reached, we 
record series of images at a constant frame rate. 
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Figure 1: Illustration of the experimental set-up: two breadboards 
mounted at 90° with the vertical one hosting most of the elements of 
the setup (SLD, camera, achromats, mirrors, polarizers, 
optomechanic items, tubing, cell, inclinometer). 


The NEFs are studied using a dynamic differential algorithm 
(DDA) that computes the spatial Fourier transform of image 
differences for each pair of images (Croccolo et al. 2012). 
Thus, we obtain the structure functions (SFs) containing 
information on the static and dynamic behaviour of the 
process. When the thermal gradient is destabilizing, the 
convective rolls emerge inside the bulk and prevent us from 
analysing the NEFs. In this case, we focus on the temporal 
evolution of the contrast of the patterns, quantified by the 
computation of the variance of the images. 


Ground-based results 

Preliminary investigations have been performed to define the 
experimental conditions of parabolic-flight experiments. It 
was found that the strength of convective rolls increases with 
increasing CO2 concentration and pressure. 

When xco, is smaller than about 10%, no convective rolls 
could be observed on ground independently of the distance 
from the bubble pressure. Nonetheless, NEFs feature a very 
limited optical contrast, so that smaller concentrations were 
discarded. For xco, = 40%, we rapidly observed convection 
for pressures in the range 6 - 24 MPa, thermal gradients of 15 
and 18 K and a mean temperature of 303 or 313K. We 
observed convective rolls that are more contrasted at high 
pressures and vanish at low pressure. These results suggest 
that at these temperatures, the Soret coefficient increases with 
decreasing xco, and approaching first-order fluid-phase 
transition. 


Results during parabola 
During the CNES parabolic flight campaign VP161, we 
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investigated two molar fractions, xco, = 26 and 56 %. We 


applied both positive and negative thermal gradients in the 
range -6 to +6K (the negative sign means heating from 
below). At the largest molar fraction, we observed convection 
for every thermodynamic condition during the terrestrial 
gravity phase. The contrast of the patterns increases during 
the hyper-gravity period, while it decreases in micro-gravity 
(Fig.4). 


Xco,= 26 % 
OT=-2K 
Trnean = 298 K 
p= 9.46 MPa 


Transition - Transition 


Xoo, = 56% 
OT=-2K 
Twnean = 297 K 
p=8.47 MPa 


Transition Transition 


Figure 2: Shadowgraph images of two experiments for a 
temperature difference AJ=-2K and two molar fractions of 
XO, = 26 (top) and 56 % (bottom). 


In some experiments convection was totally removed during 
the reduced-gravity phase; see top row of figure 2. This 
provides us periods to study the NEFs. Unfortunately, the SFs 
computed with the DDA do not provide appropriate data. 
Indeed, the vibrations of the plane and the g-jitters induce 
large perturbations of the thermally stratified layers inside the 
sample. This is especially amplified during the transition from 
hypergravity to micro-gravity. 

Despite this, we are still able to provide quantitative results 
from the analysis of the image variance. The variance of the 
raw images (var2) provides information on the presence of 
local variation of the contrast, i.e. convective patterns. The 
variance of the difference of successive images (var/) 
provides information on the movement of the patterns. 

We observe that the convective patterns are stabilized during 
the reduced gravity phase, which is illustrated with the 
decrease of var/ in figure 3. In the same period, we observed 
that the patterns are losing contrast and tend to vanish 
diffusively, which is demonstrated by the decrease of var2 in 
figure 3. From this, we are able to link the characteristic time 
decay of var2 with the thermal diffusivity of the mixture. 
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Figure 3: (left) Sample images during a parabola. (right) Temporal 
evolution of the variance of the difference of the image with the 
previous one (var/) and the variance of raw images (var2) over an 
entire parabola. The blue area corresponds to the terrestrial gravity 
phase, the green one to the hyper gravity phases, and the orange one 
to the reduced gravity phase. 
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Conclusions 

We developed an experimental set-up to study the 
establishment of convection in a mixture of COz2 and 1- 
hexanol stressed by a thermal gradient. We could highlight the 
effects of the gravity levels, from 1.8g to close to 
microgravity conditions, on the intensity of convection. This 
can vary depending on the molar fraction of the mixture and 
on the applied thermal gradient. During reduced gravity 
phases, we observed that the convective patterns are 
vanishing diffusively. The analysis of the transitory allowed 
us to measure the value of the thermal diffusivity that is in 
good agreement with literature. 

However, we were not able to study the NEFs because of the 
strength of g-jitters and vibrations of the plane. The latter 
induce a lost of the autocorrelation plane of the NEFs. In 
further campaigns, we would need to develop a better 
vibration removal apparatus, or perform the experiments in 
other microgravity facilities, like e.g. sounding rockets or 
sattelites. 
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Introduction 

The weakly nonlinear stability analysis of plane Poiseuille 
flow when viscous dissipation is taken into account is 
considered as a second part of our contribution to the IMT 
15. As we are dealing with a binary fluid mixture, the aim is 
to study the effect of thermodiffusion on the nonlinear 
properties of the instability near the onset of convection. The 
impermeable lower boundary is considered adiabatic, while 
the impermeable upper boundary is isothermal. The linear 
stability of this problem has been carried out in [Ali Amar et 
al. 2022] by solving the eigenvalue problem numerically 
using the shooting method for arbitrarily oriented convective 
rolls. The results indicated that the most unstable 
disturbances are longitudinal rolls independently of the Soret 
effect. Therefore the current study focuses on the near- 
threshold behaviour of longitudinal rolls by using a weakly 
nonlinear analysis. We determine numerically up to third 
order the coefficients of the Landau amplitude equation and 
investigate in detail the influences on _ bifurcation 
characteristics of the different dimensionless parameters 
present in the system. These parameters are namely, the 
Reynolds number Re, the Prandtl number Pr, the viscous 


dissipation Rayleigh number A defined as A = Ge Re’ Pr’ 
with Ge the Gebhart number, a parameter that measures the 
intensity of the viscous dissipation, the separation ratio y 


which represents the ratio between the mass 
contribution and the temperature contribution to 
buoyancy forces and the Lewis number Le which is 
defined here as the ratio between the thermal diffusivity and 
the solutal diffusion coefficient. In the case of a mono- 
constituent fluid (w = 0), [Requilé et al] carried out a linear 


and a weakly nonlinear stability analysis of the viscous 
dissipation induced instability both for Poiseuille and 
Couette flows. Consequently, the current study may be seen 
as an extension of [Requilé et al] to binary fluids. 


Amplitude equation of the most unstable mode 


The linear stability analysis can be used to predict when 
a particular flow becomes unstable, thus providing a 
marginal stability curve and critical values at the onset of 
the instability. However, linear stability theory cannot 
describe the evolution of unstable modes above the 
instability threshold. Nonlinear theories are used in 
order to study the further development of the instability. 
The spatiotemporal dynamics of complex flows is often 
described in fluid mechanics by relatively simple 
amplitude evolution equations. One of the most popular 
dynamic models is the Landau equation governing 
weakly nonlinear instability. This is of course not the 


only one. The attractiveness of the Landau equation is 
based on the fact that is a relatively simple model that 
allows the inclusion of physical effects. Without entering 
into mathematical details, we have proceed as in 
(Requilé et al) to derive amplitude equation that 
describes the evolution of the vertical velocity 
component in near threshold region (i.e. the viscous 
dissipation Rayleigh number A near its linear critical 
value A, determined in [Ali Amar et al. 2022]. For 
positive separation ratio y this equation reads, 

y oD De § 4 AP (1) 

ot Dee 


in which appear the characteristic time of instability v 


and the Landau coeffcient 1. Computations indicate that 
the coefficient y is real and always positive, while the 
nonlinear coefficient / is real and is positive or negative 
depending on the values of yw, Le, Pr and Ge 


considered. 
The stationary amplitude of convection reads, 


foal 
=|———+ (2) 


|A, 


Rood. 


Results of the weakly nonlinear stability 

1. Supercritical/subcritical bifurcation 

The sign of A determines whether we are dealing with a 
subcritical or a supercritical bifurcation. For A > 0, the 
bifurcation is supercritical, otherwise it becomes subcritical. 
Figure 1 shows the dependence of A on y for different 
values of Ge with Pr =10 and Le = 100. 
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Figure 1:/ asa function of y at Pr=10, Le =100 for: 
(1) Ge =10~ , 2) Ge =10", 3) Ge =1, (4) Ge=2 


As can be seen from figure 1, for reasonable values of the 
separation ratio 0<y~<0.1, the Landau coefficient Ais 


always positive for Ge = 107 indicating the supercritical 
nature of the bifurcation. Moreover, A increases with 
yw which means, according to Eq. (2), that the stationary 


amplitude of convection decreases with the increase of the 
Soret effect. However, for Ge = 107! , Ge=land Ge=2, the 


sign of A may change and the system observes a transition 
from supercritical to subcritical bifurcation at a particular 


value of y, say w=w (Le,Ge,Pr) . From figure 1 we 
observe that the more we increase the Gebhart number, the 
more the y decreases. For these prescribed values of Ge, 


figure 1 also shows that A decreases with y indicating that 


thermodiffusion reinforces the convective motions, and 
hence increases the heat and mass transfer in the fluid layer. 
As the Gebhart number plays a key role in determining the 
supercritical or subcritical nature of the bifurcation, it is 
necessary to consider its realistic numerical values. From the 
physical point of view, high values of Ge can occur only in 
geophysical or in atmospheric flows, so one can expect the 
occurrence of subcritical bifurcation for such flows. In 
contrast to these flows, if the fluid flow is on a typical 


laboratory scale, then Ge can hardly be greater than 10°. 
Results of practical relevance to experiments are presented 
in Table 1 for liquid mixtures involving ethanol/water and 
methanol/water for which the thermo-physical properties are 
well documented respectively in [Platten et al.] and 
[Lizarraga et al.]. As can be seen from Table 1, the minimum 


viscous dissipation parameter Ge = Ge(A = 0) needed to 


trigger a subcritical bifurcation is Ge =0.2316 for 


ethanol/water mixture and Ge = 0.4902 for 
methanol/water. As these values are not realistic in real 
experiments, we expect that for these particular mixtures, the 
convection may set up via a supercritical bifurcation. 


liquid properties | Water-Ethanol n-dodecane/n-hexane 


Ct 0.10 (water) 0.50 
T: (C) 37.5 25 
Pr 16.5 8.1 
Le 59.71 32.5 
u 0.071 0.076 


Nonlinear instability results 


Ag 27.49 44.36 
4 0.267 0.216 
aN 0.5788 — 2.4988 Ge 0.4753 — 0.9697 Ge 
Ge*(\ = 0) 0.2316 0.4902 
Table 1 


JI. Nonlinear effects on the convective spatial distribution 
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To illustrate the influence of the nonlinearities on the spatial 
distribution of the emerging longitudinal rolls, figure 2 
presents the streamlines in the (y, z) plane at the onset of the 


instability (Fig. 2a) and at 10% far from the criticality (Fig. 


2b). The parameters are y = 0.01, Le =100, Ge= 10“ and 
Pr =10 and the length in the y direction is fixed to the 
wave-length of longitudinal rolls at yw =0, ie. two rolls . 


As seen from Fig. 2a, the number of rolls decreases from 2 
for y =0 to 1.5 for y = 0.01. Moreover, according to linear 


theory, Fig. 2a shows a perfect mid-plane reflection 
symmetry of convection cells. In the nonlinear domain, it 
can be seen from Fig. 2b that the mid-plane symmetry is 
broken and the center of the convection rolls is shifted 
towards the upper plate, while the fluid motion is 
significantly reduced near the lower boundary. We also note 
that the nonlinearities increase the number of rolls. 


sa] 14) 


(a) (b) 
Figure 2: Streamlines: (a) Linear theory, (b) weakly nonlinear 
approach 


Conclusions 

A weakly nonlinear stability analysis of mixed convection in 
a plane parallel channel has been carried out for positive 
values of the separation ratio yw . The effect causing the 


unstable thermal stratification and, as a consequence, the 
emergence of convection cells is the viscous dissipation. 
This effect is the sole cause of the vertical temperature 
gradient in the basic state, as no thermal forcing is imposed 
through the boundary walls. We determined numerically up 
to third order the coefficients of the Landau amplitude 
equation and discussed the bifurcation nature. The main 
result obtained through this analysis is that the viscous 
dissipation, the Lewis number and thermodiffusion act in a 
concert to promote a transition from supercritical to 
subcritical bifurcation. However, for realistic values of the 
parameters in real laboratory experiments, we found that the 
bifurcation in mixed convection of binary fluids induced by 
viscous heating is always supercritical. In this case, some 
nonlinear properties of convection cells were identified in 
this work. 
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Introduction 


Water treatment technologies are at the forefront of scientific 
research to address the problem of water scarcity, which 
represents one of the worldwide unsolved issues today. 
Membrane distillation (MD) is a non-isothermal water 
treatment technology of emerging interest because of its 
environmentally-friendly benefit other than its ability to treat 
extremely high saline solutions. The MD driving force is the 
transmembrane temperature difference which induces the 
water vapor pressure difference through a _ porous 
hydrophobic membrane (Khayet and Matsuura 2011). One of 
the main challenges of MD is the design and preparation of 
membranes with suitable characteristics and improved MD 
performance (L. Garcia-Fernandez et al. 2015). Nonsolvent 
induced phase separation (NIPS) is the most commonly used 
technique for the preparation of polymeric membranes. 
However, a greener alternative is needed to mitigate the 
involved contamination of wastewaters during membrane 
preparation due to the discharge of toxic solvents, commonly 
used in industrial membrane fabrication. Membrane 
engineers are currently endeavoring to find proper green 
solvents (i.e. problematic and recommended) such as Triethyl 
phosphate (TEP), y-Valerolactone (GVL), Cyrene (CYR), 
Dimethyl sulfoxide (DMSO), and Ethyl acetate (ETAC) 
among others to replace the most toxic ones (i.e. hazardous) 
such as .N-Methyl-2-pyrrolidone (NMP), N,N- 
Dimethylacetamide (DMAC), N,N-Dimethylformamide 
(DMF); classified according to the CHEM21 guideline 
criteria (see Table 1) (Jiang et al. 2020, Wang et al. 2019, Prat 
et al. 2016). 

MD membranes must meet some exigent requirements in 
order to achieve excellent separation performance. Therefore, 
the membrane morphology should be specifically designed 
for MD, which represents a current breakthrough in the field 
of membrane science. In this research study, polymeric flat- 
sheet membranes were prepared via NIPS using the eco- 
friendly solvents shown in Table 1 together with a prior 
analysis of the thermodynamics and kinetics mechanisms. 


Table 1. Classification of solvents according to the CHEM21 
uideline criteria (Prat et al. 2016) 


Safety Health Environmental 


Solvent Ranking 


score score score 


NMP 
DMAc 
DMF 
TEP 
GVL 
CYR 
DMSO 
ETAC 


Membrane preparation via NIPS 


Firstly, an homogeneous polymer solution was prepared by 
dissolving the hydrophobic copolymer, poly(vinylidene 
fluoride-co-hexafluoropropylene) (PVDF-HFP) in different 
green solvents (DMSO, TEP, ETAC, CYR and GVL) and 
some mixed solvents. Flat-sheet membranes were prepared by 
casting followed by immersion in a coagulation bath 
composed by pure water and different nonsolvent aqueous 
mixtures (60 wt% ethanol (EtOH60), 40 wt% 
tetraethylenglycol (TetraEt40)). Finally, the membranes were 
dried at room temperature before characterization. 

The effects of the solvent type and the coagulation bath 
composition on the membrane structure were studied. 


NIPS thermodynamics and kinetics analysis 


Looking for the best MD performance (i.e. efficient heat and 
mass transfer through the membrane), a global knowledge of 
the membrane morphological structure is needed. Therefore, 
before membrane fabrication, the polymer solutions must be 
characterized by studying the polymer (P)/solvent 
(S)/nonsolvent (NS) interactions (calculated via Hansen 
solubility parameter distance, Rysp), viscosity, together with 
both NIPS thermodynamic and kinetic experiments in order 
to analyze and correlate the relationship with the membrane 
morphological structure (L. Garcia-Fernandez et al. 2014, 
2017). 

The P/S affinity is determined by Rysp (P-S), which informs 
about the solvent quality for PVDF-HFP. The solvent power 
increases for lower Rysp (P-S) values, and consequently the 
viscosity of the dope solutions decreased. This effect was 
observed when 40 wt% of TEP or ETAC was added to DMSO 
in the solvent mixture, respectively (see Fig. 1). 
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Figure 1: Rusp (P-S) and viscosity of the PVDF-HFP casting 
solutions prepared with DMSO, DMSO 60 wt%/TEP 40 wt%, and 
DMSO 60 wt%/ETAC 40 wt%. 


The NIPS thermodynamic analysis was determined by 
turbidimetric titration method. The maximum quantity of the 
nonsolvent admitted by the polymeric solution is 
denominated cloud point. It is represented in the ternary phase 
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diagram (Fig. 2) and it is a useful parameter to understand the 
thermodynamic behavior at the beginning of the phase 
inversion process occurred during membrane formation. The 
kinetic experiment showed the S-NS exchange rate during 
membrane formation. This can be studied via the evolution of 
the light transmittance during the formation of the flat-sheet 
membrane. 

By means the polymer solution characterization and NIPS 
thermodynamics and kinetics analysis, it is possible to predict 
the resultant membrane morphological structure, selecting 
therefore, the most adequate green solvent. 
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Figure 2: Ternary phase diagram of PVDF-HFP/DMSO mixed 
solvent/water system. 


Membrane characterization and water’ treatment 


application 


Various characterization techniques are used to analyze the 
suitability of the prepared membranes for MD. These are the 
thickness, contact angle, morphological structure (SEM 
cross-section and top surface), pore size distribution and 
liquid entry pressure (LEP). 

By keeping constant the solvent used for the casting solution 
preparation (DMSO 100wt%), and changing the coagulation 
bath composition, it is possible to study the coagulation 
mechanism of the different nonsolvent mixtures during the 
NIPS process. Rszp (S-NS) informs about the S-NS 
interaction, then, lower Rsxp (S-NS) values mean higher S-NS 
affinity, leading to different mean pore sizes and pore size 
distributions. As it can be seen in Fig. 3, the mixtures EtOH60 
or TetraEt40 acted as weaker nonsolvents than water (lower 
Rsup (S-NS) values). The reduced coagulation power induced 
the formation of membranes with narrower pore size 
distribution and smaller mean pore size. 
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Figure 3: Cumulative filter flow of the membranes determined by 
porometry method, and the corresponding Risp (S-NS) parameters 
obtained for different nonsolvent mixtures. 
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The morphology of the prepared membranes as well as their 
morphological structure properties can be _ explained 
according to the thermodynamics and _ kinetics 
characterization analysis. According to the membrane 
characteristics, the most suitable ones were tested in 
desalination by direct contact MD (DCMD). The obtained 
results (permeate flux and salt rejection factor) could be 
explained well with the analyzed membrane characteristics. 


Conclusions 

Polymeric flat-sheet membranes were successfully prepared 
using green solvents. These membranes exhibited suitable 
properties for MD desalination. The adopted approach, for the 
selection of the adequate green solvent through polymeric 
membrane characterization together with NIPS 
thermodynamic and kinetic analysis, proved to be interesting 
to achieve adequate MD membranes with a sustainable water 
management. 
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Introduction 


Thermodiffusion in liquid alloys influences the homogeneity 
of doped semiconductors and grown crystals as well as the 
microstructure formation during solidification. The 
calculation of thermodiffusion by molecular dynamic 
simulation can be very sensitive to the specific potential 
(Levchenko et al. 2016). Recently improved theories were 
proposed but the database to validate models for 
thermodiffusion in liquid alloys is scarce. Therefore, reliable 
measurements are needed. 


For binary alloys with sufficient contrast in x-ray absorption 
x-ray radiography can be used for in-situ measurement of 
thermodiffusion (Sondermann et al. 2019). This technique 
excludes disturbances by solidification and reveals possible 
error sources as e.g. free surfaces. 


It has been reported that in liquid Al-Cu the copper migrates 
to the cooler end of the sample (Bhat 1973) and in Al-Ag the 
silver atoms migrate to the cooler end (Kriiger et al. 2023). 
Here, we study thermodiffusion in the ternary Al-Cu-Ag melt. 


Method and Results 


In a ternary alloy it is not possible to determine the 
concentration distribution from the grey values of the x-ray 
image. It is necessary to analyze the concentration 
distribution of the processed samples after solidification. 
However, x-ray radiography is employed as process control. 


The setup for the measurement of thermodiffusion in 
multicomponent liquid alloys is very similar to the setup 
described in (Kriiger et al. 2023). It consists of a furnace made 
out of graphite and boron-nitrite which is placed inside a 
vacuum chamber. The alloy samples have the form of a 
cylinder with 1.3 mm in diameter and 10 mm in height. The 
crucible that holds the samples in vertical position is made of 
boron-nitrite. It can be separated into six parts of 1.5 mm 
height. 


During an experiment, the Al-Cu-Ag samples in the furnace 
are heated up to 973 K with a temperature gradient of 1.5 
K/mm along the sample. After annealing for at least three 
hours the liquid alloy is separated into six parts and then 
cooled down to room temperature. The solidified parts are 
subsequently removed from the crucible and analyzed by 
energy dispersive X-ray spectroscopy (EDS). 


For the eutectic composition Ales 6Agi7.6Cui3.s it is found that 
aluminum migrates to the hot end and silver migrates to the 
cold end. The copper concentration stays about constant along 
the sample. The concentration gradient of silver is 


comparable to measurements of thermodiffusion in liquid 
AlsoAgz0. 


Although the heaviest element in the alloy tends to diffuse to 
the lower part of the sample and aluminum as the lightest 
element diffuses to the upper part this is not due to gravity. In 
liquid alloys the thermal energy of an atom is larger than the 
gravitational energy by orders of magnitude (Sondermann et 
al. 2022). 


Conclusions 


With the newly developed setup, it is possible to measure the 
Soret effect in multicomponent liquid alloys. In the case of 
eutectic Al-Cu-Ag it was found that aluminum diffuses to the 
hot end while silver diffuses to the cold end. 

The Soret effect in the binary liquid Ag-Cu has not yet been 
determined due to its low contrast in x-ray absorption. From 
the results in the ternary Al-Cu-Ag we expect silver to migrate 
to the cold end also in the binary Ag-Cu melt. 
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Introduction 


Thermodiffusion in polymer solutions has revealed a number 
of surprising properties, like the molar mass independence of 
the thermodiffusion-coefficient D7 in the dilute limit or the 
blindness of the Soret-coefficient against the glass transition 
at higher polymer concentrations. Little is known, however, 
about the temperature dependence of the Soret effect in 
polymer solutions. In this work we investigate how the 
diffusion- and the Soret-coefficient change with different 
ambient temperatures, thereby covering a temperature range 
from 10°C to 50°C in steps of 5°C. The investigated system 
is a binary mixture of polystyrene in toluene. In a similar 
concentration the same system was also used in the DCMIX4 
microgravity campaign on ternary mixtures and is currently 
investigated in the ongoing GIANT FLUCTUATIONS 
project of ESA. 

The measurements have been performed with a Thermal- 
Diffusion-Forced-Rayleigh-Scattering (TDFRS) _ set-up 
(Fig.1). With this heterodyne procedure one is able to directly 
get information of the changes in the refractive index through 
the measured signal (J. Rauch et al. 2003). 


Experimental set-up 


For the TDFRS set-up we need two laser beams with different 
wavelength. One is the writing beam with a wavelength of 
532nm, which produces a holographic grating with a period 
of ~10um inside the cell. With the addition of a negligible 
amount of an inert dye, which is compatible with the 
wavelength of the laser, absorption occurs and temperature 
differences well below 0.lmK appear. This periodic 
temperature grating leads to a superposed concentration 
grating and both give rise to refractive index gratings. The 
readout laser is 632.8nm where it is important that the 
wavelength is not absorbed by the dye. The signal is 
measured by Bragg-diffraction (W. Kohler et al. 2000). 


laser A=532nm 
laser A=632.8nm 


Figure 1: The TDFRS set-up. 


Experiment 


The liquid is a binary mixture consisting of polystyrene (PS) 
with a molar mass of 4840 g/mol dissolved in toluene (Tol), 
which is a well established and researched system. The weight 
fraction of PS is cps=1%. 


We have measured the diffusion- and Soret-coefficient for 
different temperatures (10°C, 15°C, 20°C, 22°C, 25°C, 30°C, 
35°C, 40°C, 45°C, 50°C) and compared them to literature 
data at room temperature (J. Rauch et al. 2003). Both can be 
seen in Fig.2 and Fig.3 respectively. 
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Figure 2: The diffusion-coefficient of the sample PS/Tol with cps=1% for 
different ambient temperatures. The measured values of the diffusion- 
coefficient are the red dots, the literature data (J. Rauch et al. 2003) is a black 
diamond. 
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Figure 3: The Soret-coefficient of the sample PS/Tol with cps=1% for 
different ambient temperatures. The measured values of the Soret-coefficient 
are the red dots, the literature data (J. Rauch et al. 2003) is a black diamond. 
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Results 


We found a good agreement for the literature data (black 
diamonds, which were also measured with TDFRS) and our 
measured data points (red dots). The diffusion-coefficient 
rises linearly with an increase in temperature while the 
opposite can be said about the Soret-coefficient. Dr is slightly 
increasing. 


We also looked into the hydrodynamic radius (Rn) for the 
different temperatures. We calculated it with the Stokes- 
Einstein equation Rp=(kgT)/(6mND) (A. Einstein 1905). kg is 
the Boltzmann-constant, T the temperature in Kelvin, n the 
dynamic viscosity and D the diffusion-coefficient. We can 
approximate the value of the dynamic viscosity of the pure 
solvent because our sample is diluted (D. Stadelmaier et al. 
2008). The equation for the approximation is (F.J.V. Santos 
et al. 2006): 


In) Sop Oa 8.964 5.834 re 2.089 
n(j) = —5. —=— - = — 
: T T? Tr 


with the reduced variables 7 = no(T)/qo(To) and T =T/To, 
where To=298.15K and no(To)=554.2uPas. The values for the 
hydrodynamic radius follow a constant behaviour which can 
be seen in Fig.4, indicated with a horizontal black line whose 
value is Rn=1.51nm. As a result it can be said that the increase 
in temperature has no effect on the quality of the solvent and 
hence, on the expansion of these short polymer chains. 
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Figure 4: The hydrodynamic radius of the sample PS/Tol with cps=1% for 
different ambient temperatured. The values were calculated with the Stokes- 
Einstein equation.. 
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Introduction 


We report about measurements on  Poly(N- 
acryloylglycinamide) (PNAGA) in aqueous solution and 
binary and ternary mixtures of polystyrene and toluene with a 
scaled down double-pass Optical-Beam-Deflection (OBD) 
setup. For PNAGA, Diffusion-, Thermodiffusion- and Soret 
coefficients were measured with the polymer in linear form 
and as a crosslinked microgel. Furthermore, temperature 
dependent turbidity measurements were done with the linear 
PNAGA polymer. For polystyrene in toluene, results of OBD 
measurements were processed with a CONTIN routine. 


Experimental 
Experiments on diffusion were performed with a scaled down 


double-pass Optical-Beam-Deflection setup. The top view of 
the setup is sketched in Figure 1. 


Figure 1: Sketch of the top view of the scaled down double-pass 
Optical-Beam-Deflection setup. 


A laser beam transmitted through the sample gets back- 
reflected behind the OBD-cell and, thus, traverses the sample 
volume a second time. The cell with the path of the laser beam 
is sketched in Figure 2. Applying a vertical temperature 
gradient to the sample leads to thermodiffusion and the beam 
gets deflected by the resulting gradient in refractive index due 
to the gradients in temperature and concentration. A beam 
splitter guides the beam towards a camera, where the time 
dependent laser position is recorded. The signal then gives 
information about the diffusion processes. Since the laser 
beam propagates on the same path between the beam splitter 
and the cell before and after getting deflected in the sample 
and since the double-transmission geometry allows for a short 
distance between cell and camera, the whole setup could be 
built very compact. 


Cooling jaw 


Figure 2: Side view of the cell containing the sample in the Optical- 
Beam-Deflection technique. The path of the deflected laser beam is 
drawn in red. 


Results 


In case of mixtures of polystyrene/toluene, the multimodal 
time traces were evaluated with an adapted version of the 
CONTIN program for the solution of inverse problems. We 
found that diffusion coefficients extracted in this way agree 
with literature data (J. Rauch et. al. 2003). The results for 
diffusion coefficients in the measurements with samples with 
two polymers of different sizes are in agreement with the data 
for the single components. Based on these results, polymers 
with a broad molar mass distribution can now be investigated. 
Turbidity measurements with Poly(N-acryloylglycinamide) 
in water show a UCST-behavior and a hysteresis in the 
temperature of the cloud point, depending on the direction of 
the temperature ramp, which is in accordance with results 
from other experiments (J. Seuring et. al. 2012). 
Thermodiffusion- and Soret coefficients at different 
temperatures are described with an exponential Piazza 
function (S. Iacopini et. al. 2006) and are comparable for the 
linear polymer and the microgel, showing a behaviour and a 
change in sign in the region of about 35 to 40 °C. 
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Introduction 


Transport mechanisms due to temperature gradients have 
been widely studied in fields, such as energies (Mahamoudou 
et al. 2022) or the oil industry (Touzet et al. 2011). Quantiy 
describing this process is the Soret coefficient (S;), which is 
the ratio between thermodiffusion (D;) and diffusion (D) 
coefficients. For the case of a binary sample, the mass flux is 
defined as follows: 


J =—pDVc — pDpco(1 — co) VT (1) 


Where j is the mass flow, p is the density of the mixture, 
c is the concentration of the densest component in the 
mixture; Vc and VT are the concentration and temperature 
gradients respectively. When the mixture is in steady state 


> 


(J = 0) the Soret coefficient in free and porous media (Platten 
et al. 2004), S; is defined as follows: 


_ Dr Ac (2) 
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With the aim of studying this phenomenon in depth, several 
international projects have emerged, such as the Benchmark 
of Fontainebleau (Platten et al. 2003). Furthermore, there are 
various experimental techniques in free media such as the 
thermogravitational column (Blanco et al. 2008) or ODI 
(Mialdun et al. 2008) for quantifying these coefficients. In 
porous media, Soret cell (Costeséque et al. 2003) can be used 
for both positive and negative Soret coeficcients. Therefore, 
given the importance of thermodiffusion, the aim of this work 
is to set up the Soret experiment to determine its magnitude 
in terrestrial conditions in aqueous-salt mixtures. 


The article is structured as follows. First, both the device and 
the experimental procedure are briefly described. Then, the 
numerical and experimental validation with previously 
analysed samples is presented. After that, new aqueous 
samples based on LiBr salt are analysed, due to the 
importance of their transport properties in absorption chillers 
and heat pumps. Finally, conclusions are drawn from the 
study carried out. 


Experimental device and technique 


The experimental setup consists of two stainless steel cavities 
through which water circulates at a controlled temperature by 
means of two thermostatic baths LAUDA ProLine RP 855 and 
RC 6. Between the two plates, there is the gap (380x90x15 
mm) which is surrounded by PVC (Figure 1a). Both the 
stainless steel plates and the gap of constant dimensions must 


be perfectly horizontal to avoid disturbances in the tests. In 
addition, the Soret cell has 6 extraction outlets at the top and 
6 at the bottom (Figure 1b). To ensure the sealing of the 
working space during the tests, all the parts are fixed together 
with a viton gasket in between. 


Figure 1: a) Transverse cross-sectional view of the Soret cell. b) 
Whole installation of the Soret cell. 
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The experimental procedure starts by filling the gap and the 
porous medium with the mixture to ensure that the porous 
medium is completely saturated. To do so, the cell is placed 
in a vertical position and the mixture is circulated from the 
bottom upwards until the mixture comes out of the top 
bubble-free. After that, it is returned to a completely 
horizontal position and the walls are exposed to the 
corresponding temperatures. If the mixture has a positive 
Soret coefficient, the top wall is heated and the bottom one 
cooled. In contrast, if the Soret coefficient is negative, the 
other way around. This way, the denser component always 
migrates downwards. 


Numerical and experimental validation 


For the numerical validation, a 2D computational domain is 
used. The validation is carried out using a well-known 
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mixture: tetralin (THN) — isobutilbenzene (IBB) with a mass 
fraction of 0.5. The numerical model includes the 
momentum, heat and mass transfer equations in the 
Boussinesq approximation (Naumann et al. 2012). The 
density changes due to temperature and concentration are 
taken into account by the following: 


Pp = po(1 — aT — Ty) + Be — co) (3) 


Where a@ and f correspond to the coefficients of thermal 
and solutal expansion respectively. The applied temperature 
difference is 10 K. ;Error! No se encuentra el origen de la 
referencia. shows the computed concentration profile on the 
steady state. 


Figure 2: Concentration profile of THN in the THN-IBB sample 
with a temperature difference of 10 K. 
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The figure shows that the denser component (THN) goes to 
the lower part and the less dense component (IBB) 
accumulates at the upper part. 


For the experimental validation samples of 3 mL are extracted. 
In particular, the fluid from three extraction sockets is 
analysed together, obtaining four samples (two from the upper 
part and two from the lower part). In this way, enough liquid 
is obtained to calculate the concentration based on density, 
using a high-resolution density-meter. In addition, by having 
samples from two different points of the column, the values 
obtained are averaged, thus minimising experimental error. 


Following this procedure, samples previously analysed by 
(KGniger et al. 2009) have been used to perform the validation. 
The results obtained for a pair of mixtures THN-IBB and 
water-ethanol (H2O-Eth) are shown in Table , as well as the 
temperature difference and the separation of the components 
as an example. 


Table 1: Experimental results of experiments performed with the 
validation samples. 


Sample c, AT Ac S,/103 = -$,/10° 
CO) (K:') (K:!) lit 
THN-IBB 0.5 14.0 0.0123 3.49 3.3 
H20-Eth 0.9 -15.4 0.0093 -6.73 -6.75 


H20-LiBr samples 


Once the experimental device has been validated, the sample 
LiBr-H20 has been measured in the temperature range of 30 
to 45 °C and mass concentration range of 0.5 to 0.65. This 
mixture plays an important role in the field of renewable 
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energies. 

Before determining the transport coefficients of these samples, 
their thermophysical properties have been studied first. Then, 
the sliding symmetric tubes technique will be used to 
determinate their diffusion coefficients, and their Soret 
coefficients will be determined using the Soret cell with 
porous media. 


Conclusions 


The validation experiments show that the Soret cell with 
porous media is an adequated technique to measure the Soret 
coefficient of both positive and negative samples of binary 
mixtures in terrestrial conditions. 
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Introduction 


Marangoni convection greatly affects the material’s 
processing since it enhances heat and mass transport near a 
free surface. Its benefits are clearly noticed in Space due to 
the lack of natural convection but also on Earth, as seen in 
crystal growth techniques. 

Currently, the effect of thermocapillary flow in heat transport 
during a phase change material (PCM) has risen a lot of 
attention. Moreover, the use of PCMs for the extraction or 
storage of heat is of great interest, not only for terrestrial 
applications but as well in Space, mainly due to the simplicity 
of its design and usage. In this context, the literature shows 
different studies of the fluid dynamics during the melting and 
solidification phenomena, as two separate processes (see Seta 
et al. 2021 and Lappa et al. 2017) and very scarce studies that 
analyze the complete melting-solidification cycle. This could 
be interesting due to the fact that these PCM devices are 
foreseen to work continuously with this kind of cycles. Some 
of the typical materials that already have been studied during 
melting or solidification processes are sodium nitrate 
(NaNO3, Pr=8), succinonitrile (SCN, Pr=23), or the 
conventional Gallium (Ga, Pr=0.02) (Hannoun et al. 2003, 
Lappa et al. 2017). 

Paraffin-based PCMs, such as n-octadecane, are attractive 
materials due to their high latent heat values. Furthermore, 
their low density and chemical stability make them 
successfully compete with other materials already in use for 
this purpose. From this perspective, the MarPCM ESA project 
investigates the effectiveness of Marangoni convection to 
improve the heat transfer of passive devices based on PCM 
materials, Porter et al. 2022. Pure and blended n-octadecane 
under thermodiffusive Soret effect conditions are investigated 
in this project (Sanjuan et al. 2022). This work aims to study 
the flow dynamics of the n-octadecane during a complete 
melting-solidification cycle. In addition, a comparison 
between n-octadecane and the typical PCM materials 
specified earlier is proposed. To do so, 2-D simulations in a 
rectangular box were planned and the patterns of the flow 
were analyzed. 


Methodology 


The geometry used in the simulations was a rectangular cell 
with the height of 0.8 mm and the length of 8 mm (aspect ratio, 
AR=10). Initially, the material selected was in a solid state 
with a temperature, Ts, one degree under the melting one (see 
Table 1). By applying a 20 K temperature difference, AT, 
between the opposite walls (with the high temperature above 
the solid one, at the right wall) the system began to melt. The 


melting process continued until reached a minimum of 98% 
of the liquid phase. The high temperature is then lowered to 
Ts-AT, when the solidification process began. The upper wall 
was considered a free surface, where thermocapillary 
convection appeared, driving the liquid phase flow. All the 
boundary walls were assumed adiabatic and no-slip fluid 
velocity conditions were imposed except in the upper free 
border. All the simulations were carried out by using 
OpenFOAM tools. Entalphy-porosity method was used to 
model the melting-solidification cycle. In this case, the latent 
heat generated during the phase change was introduced as a 
new term in the energy equation. More details in Seta et al. 
2021. For the discretization of the system, we have checked 
different meshes and selected a 260 x 40 points mesh refined 
near the free surface and lateral walls, where high gradients 
of the physical magnitudes were expected. The physical 
properties of all the materials used in the simulations are 
summarized in Table 1. 


Table 1: Physical properties of materials used in this work: Seta et al. 
(2021)!, Lappa et al (2017) 2, Hannoun et al (2003). 


N-octad.! | SCN? NaNO,” Ga? 
Melting 
jer eo) SOLS SSIS} 581.15 302.93 
Density 5910 
tae ss 865 998 2261 
solid/liquid 6093 
Kg/m? 780 998 1908 
Specific latent heat 
(KI/Kg) 243.5 46.9 178.0 80.2 
Specific heat capacity 1934 1955 1780 37] 
solid/liquid 
Wg K) 2196 2000 1655 382 
penance ty 0.36 0.23 0.68 41 
solid/liquid ; : 5 
Wim K) UNE 0.23 0.52 32 
Dynamic viscosity 
(x 10?Pa.s) 3.6 2.6 2.7 1.8 
Thermocapillary 
coefficient 8.44 10.63 5.83 8.99 
x 10° N/(m K) 
Prandtl number 56 23 8 0.02 


Results 


Figure 1 shows the evolution of the liquid fraction along a 
melting-solidification cycle for all materials analyzed. Notice 
that the melting time of the n-octadecane is substantially 
increased (18000 s) with respect to the melting time of the rest 
of the materials (around 1800 s). This fact may be due to its 
high latent heat value which increases the melting time as well 
as its low conductivity and higher dynamic viscosity (high 
Prandtl number) that contribute to this effect. 
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In addition, the melting and solidification processes showed 
to be asymmetric, presenting higher solidification time 
compared to the melting one. During the cooling process, the 
four samples exhibited very different behavior, Gallium 
having the lowest solidification time. 
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Figure 1: Liquid fraction evolution during melting/solidification 
cycle for different materials. 


Figures 2 and 3 plot the flow pattern during the melting and 
solidification processes for all the materials, respectively. In 
both cases, the screenshots correspond to the time at which 
the liquid fraction reached 0.8 (80% melted). 


Figure 2: Screenshots of flow pattern at 0.8 liquid fraction: A) N- 
octadecane, B) SCN, C) NaNOs3, and D) Gallium. 


We can differentiate three types of flow patterns. For n- 
octadecane and SCN, the longitudinal traveling hydrothermal 
waves were formed at the hot wall (right side), with multiple 
vortex cells (see Figs. 2a and 2b). The same traveling 
hydrothermal waves were seen in the case of NaNO3 material 
though they were generated in the upper corner of the cold 
wall (left side, see Fig. 2c). In case of Gallium no traveling 
waves were detected, hence two stationary vortex cells were 
formed (see Fig. 2d). In addition, the above-mentioned flow 
patterns, due to Marangoni convection, as well conditioned 
the shape of the interphase liquid-solid. 

Once the solidification process started, two vortex cells 
appeared in all cases, originated in both right and left upper 
corners of the box, respectively (See Figure 3). These vortices 
presented a clockwise and counterclockwise rotation due to 
the existing Marangoni forces, which have opposite signs in 
the left and right corners. Remark the shape of the 
solidification front that remained practically perpendicular on 
the free surface in all cases. 
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Figure 3: Screenshots of flow pattern at 20% solidification: A) N- 
octadecane, B) SCN, C) NaNO3, and D) Gallium 


Conclusions 


Generally, for all materials investigated, melting and 
solidification processes were not symmetric presenting longer 
solidification times compared to those of melting. Moreover, 
the flow pattern and the shape of the liquid-solid interphase 
strongly depend on the Prandtl number. Under the specified 
simulation conditions, this behavior has been observed 
mainly in the melting process. In particular, n-octadecane 
showed to have a longer solidification time compared to the 
melting one, so that raises the question: do we need to wait 
for the complete cycle or it is possible to have an optimized 
heat extraction just by using a partially melting-solidification 
cycle? Therefore, more investigations are needed to answer 
this question. 
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Introduction 

The fullerene Ceo is an allotropic form of carbon where the 60 
atoms are arranged forming a hollow sphere, buckyball, with 
12 pentagonal and 20 hexagonal faces. Due to unique 
physico-chemical properties, the discovery in 1985 opened 
new routes in nanoscience and nanotechnology (Kroto et al. 
1985). Fullerenes are promising molecules for use in 
composites, energy materials, lubricants, etc., but its high 
stability, strength and low toxicity also make it a strong 
candidate for biomedical use (Bakry et al. 2007). Considered 
as potential drug delivery system, these carbon allotropes are 
active molecules for transport or can even capture molecules 
in their cage. 

Nevertheless, one of the most relevant properties that 
influence both the ease of extraction and purification, as well 
as the subsequent functionalisation and use of fullerenes, is 
their solubility (Ruoff et al., 1993). Even if fullerenes are the 
only allotropic modifications of carbon that are soluble in 
organic compounds, their solubility is in most cases low, and 
this is one of the main obstacles to progress in the study of 
their transport properties. Works in literature are almost non- 
existent; the diffusion of Coo through a semiconducting 
polymer film (polyfluorene) was analyzed using an optical 
method based on photoluminescence quenching (Fischer et al. 
2014), whereas the Soret and diffusion coefficients of Ceo 
fullerene in organic solvent o-dichlorobenzene used for 
photovoltaics were measured by the thermal-diffusion forced 
Rayleigh scattering (TDFRS) technique (Matsuura et al., 
2015). 

Motivated by the possibilities of fullerenes, the gap in the 
literature and the inherent need to determine their transport 
properties for the development of new applications, we 
present thermodiffusion (D;), molecular diffusion (D) and 
Soret coefficients (S;) of Ceo in five aromatric solvents: 
Methylnaphthalene, Toluene, and Xylene isomers oXylene, 
mXylene and pXylene. Although the structures of the 
molecules are similar, all contain cyclic atom arrangements 
and methyl groups, the thermophysical properties are 
different; That is why, in this work, we analyse the effect of 
the number of rings in the molecular structure (MN|Tol) and 
the number and position of the methyl group 
(TolloXyl|mXylijpXyl) over the mentioned transport 
properties. This work is considered as a continuation of the 
analysis of transport properties of the ternary mixture 
Coo|/THN|Tol within the framework of the DCMIX project, 
where abovementioned properties are being measured both on 
ground and on-board the International Space Station (Errarte 
et al., 2022). 


Investigated systems 

We clasify the mixtures in three blocks: diluted fullerene 
nanofluids Coot {MN|TolloXyl|mXyl|pXyl}, solvent-solvent 
binary mixtures and the ternary mixtures. 


In the case of the former, solute concentration range in mass 
composition, c; = 0.0010 — 0.0025 was analysed. Solvent 
combinations were tested at c; = 0.50. All studied ternary 
mixtures consisted of Ceo at a mass fraction of 0.002, MN at 
a concentration of 0.499 and each of the other compounds - 
Tol, oXyl, mXyl or pXy]l - at another 0.499. 

All the mixtures, represented in Figure 1 were analysed at a 
mean temperature of 25 °C. 


[ Position + number of | 
Single vs. double methyl groups 
cycle CH 
3 
CH; 
CH, 
+ ~. + ‘ 
Ke H “CH; 
CH, 
Fullerene Tol vs. MN | Tol vs. Xyl isomers | 


Figure 1: Schematic representation of molecules and analysis 
perfomed in this work. 


Experimental setups 

We tested all mixtures using the multiple-color Optical Beam 
Deflection (4-OBD) technique (Gebhardt and Kohler 2015). 
In order to determine the thermodiffusion coefficients, we 
performed supporting experiments on ternary mixtures by the 
thermogravitational column technique (TGC) (Larrafiaga et 
al. 2015). 

The former consists of determining the concentration 
variation induced by the Soret effect studying the deflection 
of a perpendicularly incident laser beam on a Soret cell. A 
mixture in thermal equilibrium is subjected to a temperature 
gradient. That condition induces components migration, 
which in turn generates a measurable refractive index gradient, 
as can be seen in Figure 2. 
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Figure 2: Normalized OBD signal from red laser (405.0 nm) for Coo, 
MN and oXy]1 based mixtures appliying a 1 °C difference. 
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Fitting analytical solutions to the OBD curves (Figure 2), the 
relaxation time t and the normalized amplitude MM (ratio 
between the value of the amplitude at the Soret plateau and 
the thermal plateau) can be determined, and, thus, the 
molecular diffusion and Soret coefficient of a binary mixture 
can be determined from D = h?/t and 


a M (=) eg (1) 
PB €9(1 = C0) \OT/) pc \OC/ pr’ 


respectively. Here, A is the cell height and cg the initial 
concentation. In the case of a ternary mixture, two 
eigenvalues of the diffusion matrix and two independent Soret 
coefficients are obtained. 

In the thermogravitational technique, a mixture is confined in 
a narrow slot between two plates with different temperatures. 
The imposed horizontal temperature gradient induces 
separation of composition due to the Soret effect in the same 
direction. At this point, a counterflow that tries to homogenize 
the system is risen. This horizontal separation also results in 
convective flow driven by buoyancy forces, which in the end 
leads to an enhanced concentration separation between the top 
and bottom ends of the column. Mixture is placed inside the 
column until the time of equilibrium, were samples are 
extracted from several intakes evenly distributed along the 
height of the column. Getting to know the concentration 
variation in function of the colum height, 0c;/dz, the 
thermodiffusion coefficients are determined by 

ass Le EPO 
TE 504 pw Oz’ (2) 


where L, corresponds to the thickness of the liquid inside 
the colum, p is the density, @ the derived thermal 
expansion coefficient, the dynamic viscosity and g the 
gravitational acceleration. Figure 3 shows the concentration 
variation along the TGC of mixture Ceo|MN\oXyl. 

0.6 ; ; r ; 


th aD coe: os. cis oe ies 
| _@ —- 8 ==e-=#--, 
0.4} ee 
a cMn = -0.067-+-0.531, R° = 0.999 
0 Coxyl = -0.068+0.466, R? = 0.999 
Fag CCyy = -0.001+0.002, R? = 0.947 
3 
0.0-—-—@ —e e eo oe 4 
0.0 0.2 0.4 0.6 0.8 


Column Height (m) 
Figure 3: Concentration variation in function of the height of the 
LTC for Coo|MN|oXyl 0.002|0.499|0.499 mixture at 25 °C. 


As can be seen, both techniques require the knowledge of 
several thermo-optical properties. The concentration effect 
over the refractive index (0n/0c) was determined measuring 
the refractive index of different samples by the Abbemat WR 
and MW refractometers. The thermal contrast factor 0n/0T 
was measured by interferometry. Density and derived thermal 
expansion coefficients were determined by the Anton Paar 
DMA 5000M density meter and viscosity by the automated 
Anton Paar AMVn micro-viscometer. 


Conclusions 

Comparing the properties among solvents differing in the 
number of solvent rings of the structure (MN|Tol), we 
observed a decrease in both diffusion and thermodiffusion 
coefficient for the double ring methylnaphthalene, but Soret 
coefficient is symilar to that of toluene. 
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Regarding the effect of methyl groups, xylene isomers show 
slightly lower diffusion and thermodiffusion coefficients than 
toluene. The order of the value of the coefficients changes 
with the shape of the organic compound, depending on the 
symmetry given by the position of the methyl group. We 
observed the same trend in the analysis of the solvent-solvent 
binary mixtures composed of MN+(TolloXyl|mXyl|pXyl). 
The eigenvalues of the diffusion matrix of ternary mixtures 
related to the Ceo are small, while the ones related to the 
liquids agree with diffusion coefficients of the corresponding 
binary subsystems MN+(ToljoXyl|mXyl|pXyl), which, again, 
verify the effect of the position of the methyl group in the 
molecules. The ill-conditioning of the solutal contrast factors 
hinder the determination of both thermodiffusion and Soret 
coefficients of the ternary systems, where in the case of the 
small coefficients of Coo, the sign is not maintained for all 
systems. Positive S; and Dj demonstrate thermophobic 
behaviour of MN and negative, hence thermophilic behaviour 
of the third, which is in good agreement with the coefficients 
of the associated binary sub-systems. 
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Introduction 


We present the thermodiffusion (Dr) coefficients of a 
polymer (polystyrene (ps) 4880 g/mol) in pure solvents 
toluene (tol) and cyclohexane (Ch) at room temperature 
(25°C), atmospheric pressure and different mass fractions of 
the polymer 2-5-10-15-20 %. For this purpose, all the 
thermophysical and optical properties of the mixtures were 
also measured in the Fluid Mechanics Laboratory of the 
University of Mondragon. The objective of this work is to use 
the thermogravitational column (TGC) procedure (convective 
technique) to determine Dr in polymeric samples and thus 
confirm the results obtained so far (tol solvent) using a non- 
convective method, thermal diffusion forced Rayleigh 
scattering (TDFRS) (J. Rauch et al. 2003). 


Methodology 


For the experimental analysis, a cartesian configuration TGC 
was considered. This consists of two parallel-plane plates in 
which the mixture to be studied is injected and a horizontal 
temperature difference is applied between them Figure 2 b). 
Regarding the method of analysis and determination of Dr, 
two main procedures can be differentiated. The first is the 
traditional method and consists of extracting different 
samples at distinct points of the TGC. After that, in the case 
of a binary mixture, the variation of the density along the TGC 
is determined and this is converted to concentration. By 
means of the slope it is possible to determine the Dr 
coefficient (Lapeira, E. et al. 2017). The second method of 
analysis is based on optical ditigital interferometry 
considering a thermogravitational microcolumn (micro-TGC). 
In this case, a Mach-Zehnder interferometer is used to 
generate interference patterns during the experiment. The 
obtained phase variation is related to the concentration 
difference (via the optical properties) along the micro-TGC 
(Lapeira, E. et al. 2018). As for the mentioned procedures, 
each one has its benefits and drawbacks. With respect to the 
traditional method, in general, the columns are larger and 
therefore the amount of sample and the test time is longer. 
However, the species separation is higher and thus the 
sensitivity of the measurements is improved. As for the 
optical method, the transient state can be studied thus 
providing additionaly the diffusion coefficient. As working 
part of the cavity is reduced it leads to the decrease of the 
separation. 


Regarding the experimental methodology followed in this 
work, first the necessary thermophysical properties (for each 
sample) such as density (p), thermal expansion coefficient (a), 
mass expansion coefficient (B) and dynamic viscosity (1) 
were determined. The Anton Paar DMA 5000 M equipment 


Figure | a) was used to determine the p of the samples. The 
same dispositive was considered to determine a and f. The 
density variation with respect to temperature, a, was 
determined according to Eqn. (1) To calculate a at 25°C the 
density of the sample was measured in the range from 24°C 
to 26°C (with step of 0.5°C). The variation of the density with 
regard to concentration B was calculated according to Eqn. (2) 
To determine B five samples of different concentrations were 
prepared. For example, in the case of the ps-tol 10 % mixture, 
samples were prepared at the following mass fractions 8-9- 
10-11-12 %. The dynamic viscosity of the mixtures 4 was 
measured using the Anton Paar AMVn microviscosimeter 
Figure | b). The optical constrast factors (dn/dc) of the 
mixtures were determined using the Anton Paar Abbemat WR 
MW refractometer Figure | c). 


a) c) 


Figure 1: Equipment used for thermophysical properties 
measurements: Anton Paar DMA 5000 densimeter a), Anton Paar 
AMVn microviscosimeter b) and Anton Paar Abbemat WR MW 
refractometer c). 
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Regarding TGC experiments, in this work, these were carried 
out both in the traditional TGC (mass fractions of 2-5 % of ps 
due to sensitivity) and micro-TGC based on the optical 
analysis (mass fractions of 10-15-20 % of ps due to the 
amount of sample needed). In the first case, extraction method, 
the large column available in our laboratory was used Figure 
2 a) which consists of the following dimension: length = 980 
mm (L,), gap = 1.02 mm (Lx) and width = 50 mm (Ly). In the 
TGC, the samples were extracted at 5 different points as 
shown in Figure 2) and the density of each extracted sample 
was determined. This procedure allows to obtain the 
concentration variation (dc/dz) along the height of the TGC 
and determine the Dr coefficient of the mixture using Eqn. (3) 
where Cy is the initial mass fraction of ps and g corresponds 
to the gravity force. 
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Figure 2: TGC experiment traditional method: used large column a) 
and a cross-section sketch b) (Lapeira, E. et al. 2017). 
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As for the optical analysis method, the micro-TGC existing in 
our laboratory based on digital interferometry was used 
(Lapeira, E. et al. 2018). This has the following dimensions: 
length = 30 mm (L,), gap = 0.509 mm (Lx) and width = 3 mm 
(Ly). Regarding the setup (Figure 3), the optical installation 
contains three lasers of different wavelengths, 473 nm, 543 
nm, 633 nm. At the starting point Figure 3 (1), three 
pneumatic cylinders are placed. In this way, the different 
lasers do not interfere with each other. After that, the filters to 
clean the beam and the collimation lenses for expansion are 
located Figure 3 (2). In the final part of the configuration, the 
Mach-Zehnder intermerometer is placed Figure 3 (3). Matlab 
R2021b software was used to analyse the interference patterns 
acquired in the camera. In this process, the first step was to 
obtain the phase variation(A@) along the column from the 
captured images. The optical phase A@ is directly related to 
the refractive index difference (An), see Eqn. (4) where (A) 
corresponds to the laser wavelength. Finally, using the 
measured optical contrast factor (0n/dc) the variation of the 
concentration was obtained Eqn. (5). The thermodiffusion 
coefficient Dr was calculated using Eqn. (3). 


Figure 3: micro-TGC experiment. Real optical installation. 
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Preliminary results 


This section presents the preliminary results obtained for the 
studied mixtures. Figure 4 a) shows the p variation along the 
height of the traditional TGC (sample ps-Tol 5 % at 25 °C) 
when the mixture reaches steady state regime. As for the 
micro-TGC experiment, Figure 4 b) represents the A@ 
(sample ps-tol 10 %) and in this case the three stages of the 
experiment can be observed. 
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Figure 4: TGC experiment (p variation) mixture ps-tol 5 % at 25°C 
a) and micro-TGC experiment A@ sample ps-tol 10 % at 25°C b). 


Conclusions 


The thermodiffusion coefficients Dy of ps in pure solvents Tol 
and Ch as a function of concetration were determined 
combining the traditional TGC technique and the micro-TGC 
method based on optical analysis. The thermophysical and 
optical properties were also measured. As for the 
experimental results obtained, the determination of the Dr 
coefficient is still going on. 
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Introduction 

Most available experimental data on thermodiffusion and 
diffusion relate to binary mixtures of a certain class of 
systems such as alkanes, glycols or organic mixtures. With 
the advent of convection-free microgravity environments, the 
focus of research has strongly shifted toward ternary 
mixtures. Still, there is a lack of studies, equally among binary 
and ternary mixtures, containing nanoparticles and, in 
particular, fullerene. The knowledge of the transport 
properties of fullerene and its compounds in various organic 
and biological solvents is important for applications. 
However, such experiments are rare. 

In this work, we focus on isothermal and non-isothermal mass 
transport of fullerene Ceo in two organic solvents measured in 
the framework of the last mission of the DCMIX project. The 
Soret effect of Coo/THN|Tol mixture at 0.0007|0.6000|0.3993 
mass fraction composition is investigated at four mean 
temperatures of 20, 25, 30 and 35 °C inside SODI, the 
Selectable Optical Diagnostic Instrument, installed in the 
MSG facility of the International Space Station (ISS) 
(Mialdun et al., 2019). Likewise, the mixture is characterized 
in three ground-laboratories by the Optical Digital 
Interferometry (ODI), the Optical Beam Deflection (OBD) 
and the Thermogravitational Columns (TGC). 

We recently published first findings about the three associated 
binary subsystems, in which thermophysical, optical and 
mass transport properties such as molecular diffusion, 
thermodiffusion and Sore coefficients were reported (Errarte 
et al. 2022). This work is a continuation to the entire 
characterization of the mixture by experiments in both earth 
laboratories and microgravity. 


Experimental setups 

The focus of the work is placed on microgravity experiments 
performed in SODI, in orbit since 2009. The instrument 
contains a two-laser Mach-Zehnder interferometer operating 
at 670 nm and 935 nm, which can sequentially access five 
cells with ternary mixtures. The experiments consist of two 
stages. First, in the thermal homogenization, the temperature 
of the system is stabilized to the mean temperature. Second, a 
thermal gradient is applied across the height of the cell almost 
instantaneously. That temperature gradient induces 
component separation, from which Soret coefficients can be 
obtained by 

1 
Spi = — ap hai: (1) 

A brief description of the data processing pipeline is given. 
SODI was designed to work by the phase stepping method, 
even if the single interferogram analysis can also be used. The 


first step consists of obtaining the linear phase of the signal 
from the intensity of the interferograms. In order to analyse 
the phase variation related to the concentration variation, a 
reference image is selected after the gradient is completely 
established in the cell and is subtracted from all subsequent 
images. The calculated unwrapped phases are converted into 
2D-refractive index maps by the expression n(y,z) = 
(A/2mL)Ad(y,z), where in this case the path the beam 
crosses through the liquid volume is of L = 10.0 mm. In 
order to fit RI data to analytic solutions that describe the 
separation in a Soret cell, the calculated 2D map is converted 
into a 1D-RI profile averaging the normalized refractive index 
map in the horizontal direction. 

Finally, the analytical solution for the heat and mass transfer 
problem describing the Soret separation is attempted by the 
quasi-binary fitting of data, successfully tested in previous 
DCMIX campaigns. Even if each fitting algorithm can lead to 
slightly different results, all of them provide both the Soret 
Spi & An;' and the effective diffusion coefficients (Derr). 
Both full path and differences methods described in Mialdun 
et al. 2018 are used to determine the Soret and effective 
diffusion coefficients and results are compared to those 
determined by the novel Optical Beam Deflection based data 
processing method reported in Sommermann et al. 2022. 
Figure | shows the steps followed for the first two fitting 
methods. 
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experiments. 
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In parallel to microgravity experiments, we are determining 
Soret and diffusion coefficients by OBD and ODI techniques 
and thermodiffusion coefficients by the TGC method. 

The former, OBD, consists of determining the concentration 
variation induced by the Soret effect by studying the 
deflection of a monochromatic laser beam incident 
perpendicularly on a diffusion cell exposed to a temperature 
gradient (Gebhardt and K6hler, 2015). Such concentration 
gradient, in turn, generates a measurable refractive index 
gradient. ODI, as SODI, uses digital interferometry to analyse 
the refractive index variation at 670 and 904 nm in a Soret 
cell (Mialdun et al. 2015). The main features of data 
processing are similar to those applied in SODI. Both OBD 
and ODI provide the diffusion coefficients and the Soret 
coefficient. In order to determine the thermodiffusion 
coefficients, we use the TGC (Larrafiaga et al, 2015). We 
insert a mixture between two vertical walls at different 
temperatures. That temperature gradient induces a horizontal 
separation. Because of this horizontal separation, a 
concentration gradient appears which in turn creates a mass 
flow in the opposite direction due to the phenomenon of 
molecular diffusion. Moreover, the effect of gravity generates 


convective flows that amplify the separation along the column, 


used for the determination of thermodiffusion coefficients. 
All these techniques require thermo-optical properties such as 
the solutal contrast factors, that turns out to be the most 
critical step during the evaluation of ternary data. This 
sensitivity leads to very asymmetrically stretched error 
ellipsoids for the Soret coefficients in the 2d-space of the 
independent concentration variables. 


Conclusions 

We present evaluation of the Coo/THN|Tol mixture analysed 
in DCMIX4 campaign performed on-board the ISS. We 
discuss different data processing methods in order to address 
their effect on the transport coefficients. Likewise, we show 
complete mixture characterization performed in three ground 
laboratories and techniques. Considering the diluted nature of 
the mixture, comparison to the corresponding binary sub- 
system THN|Tol (0.60|0.40 ge") confirms the quasi-binary 
behaviour of the system. 
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Introduction 


In a binary mixture at equilibrium there are small spontaneous 
local fluctuations in concentration, temperature and velocity 
fields around the average values of these variables (J. M. Ortiz 
de Zarate et al. 2006). For fluids at equilibrium these 
fluctuations are characterized by a white noise spectrum. In 
the presence of non-equilibrium conditions, such as those 
created by applying concentration or temperature gradients, 
the non-equilibrium fluctuations (NEFs) have a static power 
spectrum that diverges as a power law at large wave vectors q 


S(q) « ——; (A. Vailati et al. 1997). 
1+(4) 

These fluctuations are orders of magnitude larger than those 
occurring at equilibrium, and are strongly affected by the 
presence of gravity, which quenches them at wave vectors q 
smaller than a typical roll-off wave vector qr, (A. Vailati et 
al. 1998). 

Concentration NEFs are not just a small perturbation of a 
macroscopic process; in fact, they cause a net mass transfer 
that coincides with the usual Fickean diffusion. Concentration 
NEFs in the presence of a macroscopic concentration gradient 
are determined by stochastic fluctuations in the velocity field 
dv, which displace volume elements of the fluid into layers of 
different concentration. The concentration of the elements 
then relaxes to the local concentration of the fluid layer with 
a certain characteristic rate y(q). The relaxation is related to 
the wave vector by the relation: 


y(q) = Dq?- (: + (*2)') 


where D is the diffusion coefficient. The roll-off wave vector 
dro also appears in this law, indicating that gravity affects the 
relaxation of fluctuations. 

For this reason, space missions have been undertaken to 
characterize the fluctuations in more detail in the absence of 
gravity. In 2007, the GRAdient Driven FLuctuations 
Experiment (GRADFLEX), designed and carried out in 
collaboration with the European Space Agency (ESA) and 
NASA, confirmed that the theoretical models based on 
linearized hydrodynamics provide a convincing quantitative 
interpretation of nonequilibrium fluctuations for ideal 
systems (A. Vailati et al. 2011). 

However, since in nature and in technological processes NEFs 
occur in nonideal systems, where strongly nonlinear diffusion 
takes place with very intense concentration gradients, it is 
very useful to investigate such systems. The Giant 
Fluctuations and TechNES projects funded by ESA will study 
the non-equilibrium fluctuations that occur in fluids during 
diffusion processes under microgravity conditions. The focus 
of the experiments will be the investigation of transient and 
stationary processes in multicomponent fluids, also in the 
presence of intense concentration gradient. 


This is the context for the present work, in which we seek to 
provide an interpretation of results obtained by exploring 
NEFs occurring in a mixture of water and glycerol subjected 
to a strong concentration gradient. 

The experiments are performed in the isothermal free- 
diffusion configuration, where two mixtures of water and 
glycerol at different concentrations are initially brought into 
contact in the stable configuration where the denser mixture 
is at the bottom of the sample cell. To perform these 
experiments, we exploit a cell called Flowing-Junction Cell 
(F. Croccolo 2019), which can generate an interface between 
the two mixing phases. The strong dependence of the 
thermophysical property on the glycerol concentration results 
in a wide range of NEFs relaxation times at a fixed wave 
vector. 

To probe NEFs, a quantitative Shadowgraph technique was 
used (F. Croccolo et al. 2007). We investigate systematically 
diffusion processes between water-glycerol mixtures with 
glycerol concentrations ranging from 0% to 95%. 

We develop an empirical model for the time correlation 
function of the scattered light and compare it to the models 
traditionally used to characterize polydisperse colloidal 
systems, such as cumulant analysis and the Schulz 
distribution model (A. G. Mailer et al. 2015). 

We validate the effectiveness of the model with 
computational simulations and by its application to the 
experimental data. The mean square deviation between the 
model and the experimental data turns out to be better for our 
proposed model than for the others analysed. 

Furthermore, with our model we have been able to extract 
information about the stratification of the sample by 
calculating the distribution of diffusion coefficients involved 
in the mass diffusion occurring in the sample. This important 
result opens-up new scenarios in which the Shadowgraph 
technique can be employed to characterize systems in which 
it is essential to characterize stratification using a non- 
invasive method. Finally, we have characterized the impact of 
gravity for systems of this type; this step is critical in view of 
the space missions in which we are involved, where systems 
strongly stratified in concentration will be characterized 
under microgravity conditions. 
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Introduction 


We present experiments on diffusion and thermodiffusion of 
polymers in mixed solvents. So far, most works on 
thermodiffusion have dealt with binary systems or ternary 
mixtures of small molecules. Binary samples with polymers 
in solvents have been studied over both a broad concentration 
and polymer molar mass range (Rauch et al. 2003). Only a 
few very recent experiments measured polymers in a binary 
solvent (Bataller et al. 2017, Garcia-Fernandez et al. 2019). 
In this work, samples made of polystyrene, toluene and 
cyclohexane have been analysed using multiple-color optical 
beam deflection (OBD) and supporting thermogravitational 
column (TGC) experiments, based on the traditional 
extraction technique (Lapeira et al. 2017) and optical digital 
interferometry method (Lapeira et al. 2018). While binary 
mixtures are readily characterized by one diffusion and one 
thermodiffusion coefficient, the number of coefficients 
increases to four plus two for ternaries. The measured signals 
show three well separated modes that can be assigned to the 
thermal diffusivity and the two eigenvalues of the mass 
diffusion matrix. A first analysis supports the picture of an 
effective solvent whose internal dynamics is decoupled from 
the one of the polymer. 


Experimental setups 


For an OBD experiment, a temperature gradient parallel to 
gravity is suddenly applied to a horizontal sample which was 
in thermal equilibrium. A laser beam which travels through 
the middle of the cell perpendicular to the gravity direction 
gets deflected over time due to the concentration and, thus, 
refractive index gradient as can be seen in Figure 1. 


normalized OBD signal 


10? 10? 103 104 
t/s 


Figure 1: Ternary OBD signal normalized to the thermal amplitude 
together with the fit. The three indicated time constants t describe 
the diffusive behavior over time of the sample. The sample consists 
of polystyrene/toluene/cyclohexane (2.0/39.2/58.8 in mass fraction). 


In order to prevent convection, the sample height is chosen in 
a way to be far below the critical Rayleigh number. The 
temperature jump causes a fast steep change in the OBD 
signal due to the thermal diffusivity of the sample. At longer 
times, the signal evolution is driven by the Soret effect and 
can be characterized by a time constant for each independent 
concentration. In the case of a polymer dissolved in a mixed 
binary solvent the time constants are nicely separated. In the 
end the system reaches a nonequilibrium steady state, 
corresponding to a plateau of the OBD signal, from which the 
Soret coefficient can be calculated. 


The TGC technique basically consists of two parallel-plane 
plates. The sample is introduced into the cavity generated by 
the two walls and the fluid is subjected to a temperature 
gradient perpendicular to the gravitational field. During the 
experiment, due to the thermodiffusion effect, a horizontal 
concentration difference is generated in the column. This 
species separation activates the phenomenon of molecular 
diffusion, trying to homogenise the concentration gradient in 
the opposite direction. At the same time, the horizontal 
density gradient produces a convective flow along the TGC 
cavity, which increases the separation of species in the 
vertical direction. As for the method of analysis and 
determination of the thermodiffusion coefficients, two main 
procedures can be differentiated. The first, the traditional 
method, consists of extracting five different samples along the 
height of the TGC. After that, in the case of ternary mixtures, 
by determining the variation of the density and the refractive 
index, the vertical concentration gradient is obtained and 
therefore it is possible to determine the thermodiffusion 
coefficient for each species. The second analysis is based on 
optical digital interferometry using a thermogravitational 
microcolumn (micro-TGC). In this case, by means of a Mach- 
Zehnder interferometer, interference patterns are acquired for 
each step of the micro-TGC experiment: temperature 
homogenisation, temperature gradient application and steady 
state regime. Figure 2 shows the phase variation in a micro- 
TGC experiment for a stable (blue) and unstable (red) sample 
showing each stage of the experiment. 
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Figure 2: Phase variation over time in a micro-TGC experiment for 
a stable (blue) and an unstable (red) mixture. 
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Investigated system 


As stated above, the system consisting of one polymer and 
one solvent has been analysed for various components over 
the last decades. In order to build upon the existing research 
and keeping the measurement time decent, polystyrene (My = 
4880 g/mol) was used in combinations with toluene and 
cyclohexane over a broad concentration range, which can be 
seen in Figure 3. The polymer concentrations and solvent 
ratios are always based on the masses of the constituents. It 
turns out that the solvent ratio has a big influence on the fast 
mode and the polymer concentration on the slow mode (t2 and 
3 in figure 1, respectively). With increasing toluene content, 
the amplitude of the fast mode decreases, while for increasing 
polymer concentration at constant solvent ratio, the amplitude 
of the slow mode also increases. For the measurement series 
with a solvent ratio toluene:cyclohexane of 80:20, the 
evaluation of the OBD signal gets more difficult because the 
fast mode is hard to distinguish from the noise of the thermal 
amplitude; the ternary signal almost resembles a binary one. 
This problem was the main reason to look for other 
measurement techniques in order to provide some reference 
data. The TGC experiments proofed to be reliable, and the 
evaluation is not so much affected by noise as in the case of 
OBD. As a downside, the Soret coefficient of the densest 
component must be positive in order to be measurable. 
Preliminary experiments with a smaller micro-TGC showed 
that not all samples are stable (see red and grey squares in 
(Figure 3), so the focus was on samples with a sufficient large 
polymer concentration (blue squares in Figure 3). These 
samples were analysed at a mean temperature of 25°C while 
the faster OBD measurements could also be investigated at 
20°C, 30°C and 35°C. 


PS (My=4880Da) in toluene:cyclohexane 
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Figure 3: Measured samples in the OBD setup (open squares) and 
additionally with the TGC setup (filled squares). The y-axis indicates 
the solvent mass ratio in terms of toluene:cyclohexane. Due to 
instabilities not all samples are measurable in the TGC. 


Data analysis 


While the data evaluation is still going on, a few tendencies 
can already be seen. The fast mode, which corresponds to the 
smaller eigenvalue of the mass diffusion matrix, matches 
pretty well in the limit of small polymer concentrations with 
the diffusion constant in the pure solvent samples. The larger 
eigenvalue, on the other hand, scales in first order 
approximation linear with the polymer concentration and is 
confined in an area which is defined by the eigenvalue of the 
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binary system of polymer in the two mentioned solvents. To 
get the thermodiffusion coefficients from an OBD 
measurement, the inversion of a so called contrast factor 
matrix is needed, which is often ill conditioned and results in 
a large error propagation. The analysis of the TGC uses a way 
better conditioned contrast factor matrix leading to 
significantly less noise in the results. A definite interpretation 
of the comparison of the thermodiffusion coefficients from 
both setups is still under debate. 


Conclusions 


The OBD technique was applied to investigate ternary 
samples consisting of a small polystyrene molecule in toluene 
and cyclohexane. It was possible to measure the system over 
a broad range, starting from the binary samples to ternaries 
with a polymer concentration up to 10%. Due to problems 
during the signal evaluation, experiments were also 
performed using TGC in order to create reference data. 
Concerning the well separated diffusion coefficients, they can 
be attributed to a solvent-solvent mode and a polymer motion 
in an effective solvent. The interpretation of the 
thermodiffusion coefficient is still going on, mainly hindered 
by error propagation. 
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Introduction 


The paper is devoted to the investigation of the 
Soret-induced convection of ternary fluid in horizontal 
porous layer heated from below. The study of convection in 
horizontal porous layers is important for understanding the 
processes of heat and mass transfer in oil-bearing fields. Oil 
is a complex mixture of a large number of different fractions 
of hydrocarbons. In oil-bearing fields it is subjected to 
geothermal temperature gradient. The knowledge of 
peculiarities of heat and mass transfer in these conditions 
could allow to make the oil extraction more efficient. The 
problems of convection of multi-component fluids began to 
be considered quite recently and are not enough understood. 


Problem formulation, governing equations, boundary 
conditions and system parameters 


An infinite horizontal layer of a porous medium 
saturated with a multicomponent fluid is considered. As a 
fluid, a three-component mixture of hydrocarbons is 
considered: dodecane (component Nel), isobutylbenzene 
(component Noe2) and tetralin (component Ne3) with equal 
mass fractions of components at an average temperature of 
25°C (LI. Ryzhkov, 2013). These components are typical for 
mixtures found in natural hydrocarbon deposits. We choose 
the heaviest component (tetralin) as the solvent. This 
mixture was studied in experiments under microgravity 
conditions within the framework of the DCMIX 1 project (V. 
Sechenyh et al. 2016). The system is subjected to the gravity 
field. It is assumed that the boundaries of the layer are rigid, 
perfectly conductive and impermeable for solute. Different 
constant temperatures are maintained at the boundaries. It is 
assumed that the density of the fluid depends linearly on the 
temperature and concentrations of the components. 

The equations of the Soret-induced convection of 
multicomponent fluid in a porous medium (N.I. Lobov et al. 
2007) in the dimensionless form are: 


0 =-VP—V+Ra,(T + pC)je, (1) 
o+V-VT =VT (2) 

eX +V-Ve = Le*(v°C -1-V°7) (3) 

V:-V=0 (4) 

where Ra, = a is an analog of the Rayleigh number 


for a porous medium (hereinafter, we call it the “Rayleigh 
number”); = —Co(1 — Co) 87'B(D)“1Dz is the vector 


(pc) ¢ 
_ ny (pc)* 
Le = a(D)"* is the diagonal matrix of the Lewis numbers, 


— a < & = 
, Les, =—; ey, is the unit vertical 
D22 ™ 


of separation ratio; € = €* is the normalized porosity; 


—— a 
where Le,, = ai 
11 


* 


x 
vector. Here a = 
(pc) 


V is the filtration velocity, P is the pressure, g is the 
gravity acceleration, K is the permeability of the porous 
matrix, v is the kinematic viscosity, A“ is the effective 
thermal conductivity of the porous medium, (pc)* is the 
effective heat capacity of the porous medium, (pc), is the 
heat capacity of the mixture, €* is the porosity, D is the 
matrix of molecular diffusion coefficients, J = (1,...,1) is 
the column vector, IJ = (1,...,1) is the row vector, C = 
(Cy, ...,Cy-1) 1s the concentration vector, Dy is the vector 
of thermal diffusion coefficients. 

At the layer boundaries maintained at constant 
different temperatures the impermeability condition and the 
condition of the absence of the solute flux are imposed. 
These conditions in the dimensionless form are: 


is the effective thermal diffusivity, 


The problem was studied for fixed values of the 
parameters wp, Le,, Lez, € i, = 0.120, fe, = 
147, Le, =91, € =0.15. The separation ratio p, was 
varied in the interval [—0.150, O]. 


Linear stability of conductive state 


To solve the linear problem of the stability of the 
conductive state, the Runge-Kutta method of the 5th order 
and the subsequent search for the zero value of the 
determinant of the matrix of the fundamental system of 
solutions using the two-dimensional secant method were 
used. 

Fig. 1 shows the linear stability map obtained in the 
calculations. For heating from below only the longwave 
monotonic instability mode exists in the range yw, € 
[—-0.047;0] (curve 1). At pw, = —0.047, the longwave 
oscillatory instability mode branches off from the longwave 
monotonic instability mode and becomes more dangerous 
(curve 3). At w,=-—0.065 , the finite-wavelength 
monotonic instability mode (curve 2) branches off from the 
longwave monotonic instability mode. At w, = —0.101, the 
longwave oscillatory instability transforms into 
finite-wavelength oscillatory instability (curve 4). At py, = 
—0.113, this curve changes the slope. 
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Figure 1: Linear stability map (instability curves: 1 - monotonic 
longwave, 2 - monotonic finite-wavelength, 3 - oscillatory 
longwave, 4 — oscillatory finite-wavelength) 


Note that the onset of convection in the system 
under consideration at similar parameter values was studied 
earlier in (D. Mutschler et al. 2020), but in this work the 
longwave oscillatory mode of instability was not found. 


Nonlinear regimes of convection 


Nonlinear calculations were carried out for 
two-dimensional horizontally elongated rectangular cavity 
with aspect ratio 5xl. At all the boundaries, the 
impermeability condition and the absence of the solute flux 
were set. The constant different temperatures were imposed 
at the horizontal boundaries, and the vertical boundaries 
were assumed to be adiabatic. The stream function was 


: : : : aw 

introduced, which was determined by the relations V, = rr 
aw 

Vy =o 


The problem was solved in an unsteady formulation 
by the finite difference method. The Poisson equation for the 
stream function was solved by’ the — successive 
over-relaxation method. A uniform grid was used. The 
number of grid cells was 20 vertically and 100 horizontally. 

Modeling of nonlinear convection regimes was 
carried out for the values y,, which are characteristic points 
on the stability map: 


YW, = {-0.055, —0.080, —0.110, — 0.130} 


The calculations have shown that at ww, = 
—0.055 after stability loss by the conductive state at Ra ~ 
3.26, the regime of stationary oscillations of the longwave 
type is excited, via direct bifurcation. With the increase in 
the Rayleigh number to Ra=3.40, this regime is 
transformed into the monotonic longwave regime. If, being 
on the curve of the monotonic regime, the Rayleigh number 
is lowered, then at Ra ~ 3.10 the convection starts to 
decay. Thus, the transition between different convection 
regimes is accompanied by hysteresis. 

At yw, =—0.080, after stability loss by the 
conductive state at Ra ~ 5.26, the longwave oscillatory 
regime is excited via direct bifurcation. When the Rayleigh 


T : T T i. 
-0.06 -0.04  -0.02 
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number increases to Ra ~ 5.95, this regime passes into a 
monotonous longwave regime. If, while on the curve of the 
monotonic regime, the Rayleigh number is lowered, then a 
reverse transition to the oscillatory regime occurs. 

At w, =—0.110, after stability loss by the 
conductive state at Ra ~ 18.58, the regime of stationary 
oscillations of the three-vortex structure is excited, which, 
with a further increase in the Rayleigh number to Ra ~ 
20.00, transforms into a monotonous five-vortex regime. If, 
being on the curve of the five-vortex regime, the Rayleigh 
number is lowered, then at Ra ~ 15.50 there is a transition 
to the four-vortex monotonic regime, which exists up to 
Ra = 13.20. With a further decrease in Ra, the convection 
starts to decay. 

At W, = —0.130, after stability loss by the 
conductive state at Ra ~ 36.85, a six-vortex stationary 
regime is excited. If, being on this curve, we reduce the 
Rayleigh number, then at Ra ~ 36.00 there is a transition 
to the five-vortex stationary regime, which exists up to 
Ra ~ 34.08. With a further decrease in Ra, the convection 
starts to decay. 


Conclusions 


In the investigated range of parameters, four types 
of instability were found in the system: long-wave 
monotonic, cellular monotonic, long-wave oscillatory and 
cellular oscillatory. 

Within the framework of linear stability analysis, 
a new longwave oscillatory mode of instability was found, 
the existence of which was confirmed by nonlinear 
calculations using the finite difference method. 

The convection threshold found in non-linear 
calculations agrees with the threshold predicted by linear 
theory. 

At wy, = -0.055 , Ww, =-0.080 and y, = 
—0.110, the direct bifurcation takes place and at wy, = 
—0.130, convection is excited in a finite-amplitude way. 

Nonlinear calculations have shown that the 
transitions between different convection regimes are 
accompanied by hysteresis. 
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Introduction 

One of the variants of the structure of an oil-bearing field is a 
horizontal reservoir. The simplest model of such a reservoir 
can be a layer of porous medium (sands, limestones) saturated 
with oil and limited by hard-permeable rocks. Such a deposit 
is called lithologically limited. However, most often the 
structure of an oil-bearing formation is some rather complex 
tectonic form, where hydrocarbon-saturated layers with 
different porosity and permeability can alternate. 

Oil is a multicomponent mutually soluble mixture of 
hydrocarbons of various structures with a small admixture of 
some other substances. The complex composition of the 
liquid filling the porous medium is one of the complicating 
factors in modeling. The properties of the mixture may be 
such that a positive or negative thermal diffusion effect may 
occur under the action of a geothermal temperature gradient. 
In this case, instability can occur, leading to the appearance 
of convection. 

The aim of this work is to simulate the convection of a binary 
mixture of liquid hydrocarbons in a two-layer porous 
medium, the layers of which have different permeabilities and 
are separated either by a flat horizontal boundary or by a 
boundary that simulates a synclinal geological fold. 


Problem configuration and governing equations 

The considered mixture of dodecane (Ci2H2.) and tetralin 
(CioHi2), taken in equal proportions [Gebhardt et al. 2013], 
fills a horizontally elongated rectangular region of a porous 
medium with the dimensions L=200 m and H=100 m 
and rigid external boundaries impermeable to the matter. The 
computational domain consists of two horizontal layers with 
different constant values of permeability K, and K,, but the 


same other properties. We study two options for the shape of 
the boundary between the layers: a flat boundary that divides 
the area into two layers of equal height (Fig. la), and a 
boundary in the form of an arc of a circle with a radius of 
187.5 m, directed downwards, which imitates the boundary of 
a synclinal geological fold (Fig. 1b). 

In the simulation, we assume that the viscosity and transport 
coefficients are constant and neglect the barodiffusion and the 
Dufour effect. Calculations of the Soret-induced convection 
of a binary mixture in a porous medium can be carried out 
within the framework of non-stationary nonlinear Darcy— 
Boussinesq equations [Nield 2013]: 


1 
are ~V,-g(B; (7, -%)+Bc(C,-G)), 
0 i 
(pc) <i +(pe),¥, VI, =hV'T,, 
- OC, 2 2 
& SE 4V,-VC, = DVC, +6, (I-G) DVT, 


V-V,=0. 


Here V is the filtration rate; p is the pressure; T is the 
temperature; C is the concentration of tetralin; £4, is the 
f- is coefficient of the 
concentration dependence of the density; 7, and C, are 
is the 
acceleration of gravity; ¢ is the time; K, and K, are the 


thermal expansion coefficient; 
the initial temperature and concentration; g 


permeability of the upper and lower layers; v is the 


kinematic viscosity of the mixture; 2” is the effective 
thermal conductivity of the porous medium; (pc) is the 


effective heat capacity of the porous medium; (pc) r is the 


heat capacity of the mixture; ¢° is the porosity; D is the 
molecular diffusion coefficient of the mixture; D, is the 
i=1,2 


thermal diffusion coefficient. Indexes indicate 


correspondence to the top or bottom layer. 


H/2 


L x 


a 
Figure 1: Problem geometry. a — system of layers of equal height, 
b — system imitating a synclinal fold. 


At all external boundaries we set the conditions of 
impermeability and the absence of a diffusion flux of matter, 
at the upper and lower boundaries we impose constant 
different temperatures, on vertical ones — the absence of a heat 
flux; at the boundary between the layers — the continuity of 
pressure, heat flux and diffusion flux of matter. 

At the initial moment of time, we assume the liquid to be 
motionless, the solute concentration in the region to be 
homogeneous, and the temperature to be linearly dependent 
on the vertical coordinate. The temperature difference at the 
external horizontal boundaries corresponds to the average 
value of the geothermal gradient: 3-10 K/m. 

The calculations were carried out using the ANSYS Fluent 
software package, which implements the finite volume method. 
A spatially uniform grid with square cells was used, the grid 
spacing was 2 m. To obtain a discrete analog of the equations, 
scheme of the second order of approximation in time and the 
third in space was used. 


Numerical results 
We considered the cases where the upper or lower layer was 
more permeable. 
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An analysis of the time evolution of the maximum value of 
the stream function in the computational domain has shown 
that for any considered configuration of the two-layer system, 
a certain non-convective period is first observed (the stream 
function is practically zero), the duration of which decreases 
with increasing permeability of any layer. The emergence of 
convection is accompanied by a sharp jump in the intensity of 
motion. Then, through a series of rearrangements, which are 
accompanied by smaller jumps in flow intensity, the steady 
state is reached. As an example, Figure 2 shows the temporal 
evolution of the maximum value of the stream function in a 
two-layer system with a curved boundary between the layers 
for the cases y=K,/K, <1 (the lower layer is more 


permeable). 


No 


x» 10°, kg/(ms) 


= 0.210" 


y T : T 

0 2 4 6 8 10 
i, 104, year 

Figure 2: Temporal evolution of the maximum value of the stream 

function in a two-layer system with a curved boundary between the 

layers, y<1; gray dashed curves correspond to the case of a flat 


layer boundary and the same value of y. 


During the non-convective period, a process of separation of 
the mixture close to pure diffusion is observed: the 
concentration difference between the upper and lower 
boundaries tends to a value corresponding to the maximum 
separation of the mixture. For higher values of the medium 
permeability, the rapidly occurring convective motion 
strongly mixes the mixture, as a result of which the maximum 
separation of the mixture is substantially reduced. 

With significantly different permeabilities of the layers (the 
permeability of one of the layers is ten or a hundred times 
higher than that in the other), a “local” occurrence and 
existence of convection is observed in the cavity, i.e. the flow 
arises and develops in a layer with a higher permeability 
(Fig. 3a). With a decrease in the difference in the 
permeabilities of the layers, a weak penetration of the flow 
into the less permeable layer is observed. At close values of 
the permeability of the layers, the emergence and existence of 
convection is of a “large-scale” nature: the flow arises and 
develops throughout the volume of the system, but is shifted 
to a layer of a higher permeability (Fig. 3b). It is found that at 
close permeabilities of the layers, one-, two-, or four-vortex 
stationary flow is established. For strongly different 
permeabilities of the layers, the appearance of quasi-periodic 
or irregular oscillations of complex shape is typical; the flow 
structure in this case is one- or two-vortex. 
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Figure 3: Fields of the stream function at ¢=20000 years. a — 
y=10°,6-— y=0.2-10". 


The distribution of the solute concentration during the “local” 
development of convection and significantly different 
permeabilities of the layers is almost uniform in the layer 
where convection is observed (a layer of a_ higher 
permeability) and close to a linear distribution along the 
vertical in a layer where the flow practically does not 
penetrate (a layer with a lower permeability). A decrease in 
the permeability difference leads to a greater penetration of 
the flow into a less permeable layer. In this case, the 
concentration isolines have a characteristic finger-like shape. 


Conclusions 

The simulation performed has shown a very slow 
development of convection (at times of the order of thousands 
of years) and the formation of steady-state regimes (at times 
of the order of tens of thousands of years). It is found that at 
close permeabilities of the layers, one-, two-, or four-vortex 
stationary flow is established. For strongly different 
permeabilities of the layers, the appearance of quasi-periodic 
or irregular oscillations of complex shape is typical; the flow 
structure in this case is one- or two-vortex. 

It is shown that the "local" development of convection (the 
flow arises and exists in a layer of a higher permeability) is 
typical for layers that differ significantly in permeability (the 
permeability of one of the layers is ten or one hundred times 
higher than that in the other). With an increase in the 
permeability of the second layer, a slight penetration of the 
flow into this layer is observed. At even closer values of the 
permeabilities of the layers (less than ten times), the onset and 
development of convection is of a “large-scale” nature: the 
flow occurs in the entire volume of the system, but is shifted 
to a layer with a higher permeability. 

The emerging convective flows have low filtration rates. 
However, these low velocities are sufficient to mix the fluid 
even in the case of weak flow penetration into the less 
permeable layer and to prevent its separation. 
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